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 develop

Browse Source
This commit is contained in:
Germain Clavier
2024-12-19 13:14:38 +01:00
parent 5bec534faf aeb2190582
commit 1e968d9cb0
371 changed files with 50393 additions and 17270 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
.github/CODEOWNERS vendored
View File

@ -101,7 +101,8 @@ src/group.* @sjplimp
src/improper.* @sjplimp
src/info.* @akohlmey
src/kspace.* @sjplimp
src/lmptyp.h @sjplimp
src/lmptype.h @sjplimp
src/label_map.* @jrgissing @akohlmey
src/library.* @sjplimp @akohlmey
src/main.cpp @sjplimp
src/min_*.* @sjplimp

108
.github/release_steps.md vendored Normal file
View File

@ -0,0 +1,108 @@
# LAMMPS Release Steps
The following notes chronicle the current steps for preparing and publishing LAMMPS releases. For
definitions of LAMMPS versions and releases mean, please refer to [the corresponding section in the
LAMMPS manual](https://docs.lammps.org/Manual_version.html).
## LAMMPS Feature Release
A LAMMPS feature release is currently prepared after about 500 to 750 commits to the 'develop'
branch or after a period of four weeks up to two months. This is not a fixed rule, though, since
external circumstances can cause delays in preparing a release, or pull requests that are desired to
be merged for the release are not yet completed.
### Preparing a 'next\_release' branch
Create a 'next\_release' branch off 'develop' and make the following changes:
- set the LAMMPS\_VERSION define to the planned release date in src/version.h in the format
"D Mmm YYYY" or "DD Mmm YYYY"
- remove the LAMMPS\_UPDATE define in src/version.h
- update the release date in doc/lammps.1
- update all TBD arguments for ..versionadded::, ..versionchanged:: ..deprecated:: to the
planned release date in the format "DMmmYYYY" or "DDMmmYYYY"
- check release notes for merged new features and check if ..versionadded:: or ..versionchanged::
are missing and need to be added
Submit this pull request, rebase if needed. This is the last pull request merged for the release
and should not contain any other changes. (Exceptions: this document, last minute trivial(!) changes).
This PR shall not be merged before **all** pending tests have completed and cleared. If needed, a
bugfix pull request should be created and merged to clear all tests.
### Create release on GitHub
When all pending pull requests for the release are merged and have cleared testing, the
'next\_release' branch is merged into 'develop'.
Check out 'develop' locally, pull the latest changes, merge them into 'release', apply a suitable
release tag (for historical reasons the tag starts with "patch_" followed by the date, and finally
push everything back to GitHub. Example:
```
git checkout develop
git pull
git checkout release
git pull
git merge --ff-only develop
git tag -s -m "LAMMPS feature release 19 November 2024" patch_19Nov2024
git push git@github.com:lammps/lammps.git --tags develop release
```
Go to https://github.com/lammps/lammps/releases and create a new (draft) release page or check the
existing draft for any necessary changes from pull requests that were merged but are not listed.
Then select the applied tag for the release in the "Choose a tag" dropdown list. Go to the bottom of
the list and select the "Set as pre-release" checkbox. The "Set as the latest release" button is
reserved for stable releases and updates to them.
If everything is in order, you can click on the "Publish release" button. Otherwise, click on "Save
draft" and finish pending tasks until you can return to edit the release page and publish it.
### Update download website, prepare pre-compiled packages, update packages to GitHub
Publishing the release on GitHub will trigger the Temple Jenkins cluster to update
the https://docs.lammps.org/ website with the documentation for the new feature release
and it will create a tarball for download (which contains the translated manual).
Build a fully static LAMMPS installation using a musl-libc cross-compiler, install into a
lammps-static folder, and create a tarball called lammps-linux-x86_64-19Nov2024.tar.gz (or using a
corresponding date with a future release) from the lammps-static folder and upload it to GitHub
with:
```
gh release upload patch_19Nov2024 ~/Downloads/lammps-linux-x86_64-19Nov2024.tar.gz
```
## LAMMPS Stable Release
A LAMMPS stable release is prepared about once per year in the months July, August, or September.
One (or two, if needed) feature releases before the stable release shall contain only bug fixes
or minor feature updates in optional packages. Also substantial changes to the core of the code
shall be applied rather toward the beginning of a development cycle between two stable releases
than toward the end. The intention is to stablilize significant change to the core and have
outside users and developers try them out during the development cycle; the sooner the changes
are included, the better chances for spotting peripheral bugs and issues.
### Prerequesites
Before making a stable release all remaining backported bugfixes shall be released as a (final)
stable update release (see below).
A LAMMPS stable release process starts like a feature release (see above), only that this feature
release is called a "Stable Release Candidate" and no assets are uploaded to GitHub.
### Synchronize 'maintenance' branch with 'release'
The state of the 'release' branch is then transferred to the 'maintenance' branch (which will
have diverged significantly from 'release' due to the selectively backported bug fixes).
### Fast-forward merge of 'maintenance' into 'stable' and apply tag
At this point it should be possible to do a fast-forward merge of 'maintenance' to 'stable'
and then apply the stable\_DMmmYYYY tag.
### Push branches and tags
## LAMMPS Stable Update Release

126
.github/workflows/kokkos-regression.yaml vendored Normal file
View File

@ -0,0 +1,126 @@
# GitHub action to build LAMMPS on Linux and run selected regression tests
name: "Kokkos OpenMP Regression Test"
on:
push:
branches:
- develop
workflow_dispatch:
jobs:
build:
name: Build LAMMPS with Kokkos OpenMP
# restrict to official LAMMPS repository
if: ${{ github.repository == 'lammps/lammps' }}
runs-on: ubuntu-latest
env:
CCACHE_DIR: ${{ github.workspace }}/.ccache
strategy:
max-parallel: 6
matrix:
idx: [ 'pair-0', 'pair-1', 'fix-0', 'fix-1', 'compute', 'misc' ]
steps:
- name: Checkout repository
uses: actions/checkout@v4
with:
fetch-depth: 2
show-progress: false
- name: Install extra packages
run: |
sudo apt-get update
sudo apt-get install -y ccache ninja-build libeigen3-dev \
libcurl4-openssl-dev python3-dev \
mpi-default-bin mpi-default-dev
- name: Create Build Environment
run: mkdir build
- name: Set up ccache
uses: actions/cache@v4
with:
path: ${{ env.CCACHE_DIR }}
key: linux-kokkos-ccache-${{ github.sha }}
restore-keys: linux-kokkos-ccache-
- name: Building LAMMPS via CMake
shell: bash
run: |
ccache -z
python3 -m venv linuxenv
source linuxenv/bin/activate
python3 -m pip install --upgrade pip
python3 -m pip install numpy pyyaml junit_xml
cmake -S cmake -B build \
-C cmake/presets/gcc.cmake \
-C cmake/presets/basic.cmake \
-C cmake/presets/kokkos-openmp.cmake \
-D CMAKE_CXX_COMPILER_LAUNCHER=ccache \
-D CMAKE_C_COMPILER_LAUNCHER=ccache \
-D BUILD_SHARED_LIBS=off \
-D DOWNLOAD_POTENTIALS=off \
-D PKG_AMOEBA=on \
-D PKG_ASPHERE=on \
-D PKG_BROWNIAN=on \
-D PKG_CLASS2=on \
-D PKG_COLLOID=on \
-D PKG_CORESHELL=on \
-D PKG_DIPOLE=on \
-D PKG_DPD-BASIC=on \
-D PKG_EXTRA-COMPUTE=on \
-D PKG_EXTRA-FIX=on \
-D PKG_EXTRA-MOLECULE=on \
-D PKG_EXTRA-PAIR=on \
-D PKG_GRANULAR=on \
-D PKG_LEPTON=on \
-D PKG_MC=on \
-D PKG_MEAM=on \
-D PKG_POEMS=on \
-D PKG_PYTHON=on \
-D PKG_QEQ=on \
-D PKG_REAXFF=on \
-D PKG_REPLICA=on \
-D PKG_SRD=on \
-D PKG_SPH=on \
-D PKG_VORONOI=on \
-G Ninja
cmake --build build
ccache -s
- name: Run Regression Tests for Selected Examples
shell: bash
run: |
source linuxenv/bin/activate
python3 tools/regression-tests/get_kokkos_input.py \
--examples-top-level=examples --batch-size=50 \
--filter-out="balance;fire;gcmc;granregion;hyper;mc;mdi;mliap;neb;pace;prd;pour;python;rigid;snap;streitz;shear;ttm"
export OMP_PROC_BIND=false
python3 tools/regression-tests/run_tests.py \
--lmp-bin=build/lmp \
--config-file=tools/regression-tests/config_kokkos_openmp.yaml \
--list-input=input-list-${{ matrix.idx }}-kk.txt \
--output-file=output-${{ matrix.idx }}.xml \
--progress-file=progress-${{ matrix.idx }}.yaml \
--log-file=run-${{ matrix.idx }}.log \
--quick-max=100
tar -cvf kokkos-regression-test-${{ matrix.idx }}.tar run-${{ matrix.idx }}.log progress-${{ matrix.idx }}.yaml output-${{ matrix.idx }}.xml
- name: Upload artifacts
uses: actions/upload-artifact@v4
with:
name: kokkos-regression-test-artifact-${{ matrix.idx }}
path: kokkos-regression-test-${{ matrix.idx }}.tar
merge:
runs-on: ubuntu-latest
needs: build
steps:
- name: Merge Artifacts
uses: actions/upload-artifact/merge@v4
with:
name: merged-kokkos-regresssion-artifact
pattern: kokkos-regression-test-artifact-*

11
cmake/CMakeLists.txt
View File

@ -141,7 +141,7 @@ endif()
# silence nvcc warnings
if((PKG_KOKKOS) AND (Kokkos_ENABLE_CUDA) AND NOT (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))
set(CMAKE_TUNE_DEFAULT "${CMAKE_TUNE_DEFAULT} -Xcudafe --diag_suppress=unrecognized_pragma -Xcudafe --diag_suppress=128")
set(CMAKE_TUNE_DEFAULT "${CMAKE_TUNE_DEFAULT}" "-Xcudafe --diag_suppress=unrecognized_pragma,--diag_suppress=128")
endif()
# we require C++11 without extensions. Kokkos requires at least C++17 (currently)
@ -588,13 +588,8 @@ endif()
set(CMAKE_TUNE_FLAGS "${CMAKE_TUNE_DEFAULT}" CACHE STRING "Compiler and machine specific optimization flags (compilation only)")
separate_arguments(CMAKE_TUNE_FLAGS)
foreach(_FLAG ${CMAKE_TUNE_FLAGS})
target_compile_options(lammps PRIVATE ${_FLAG})
# skip these flags when linking the main executable
if(NOT (("${_FLAG}" STREQUAL "-Xcudafe") OR (("${_FLAG}" STREQUAL "--diag_suppress=unrecognized_pragma"))))
target_compile_options(lmp PRIVATE ${_FLAG})
endif()
endforeach()
target_compile_options(lammps PRIVATE ${CMAKE_TUNE_FLAGS})
target_compile_options(lmp PRIVATE ${CMAKE_TUNE_FLAGS})
########################################################################
# Basic system tests (standard libraries, headers, functions, types) #
########################################################################

6
cmake/Modules/FindVORO.cmake
View File

@ -21,9 +21,9 @@ if(VORO_FOUND)
set(VORO_LIBRARIES ${VORO_LIBRARY})
set(VORO_INCLUDE_DIRS ${VORO_INCLUDE_DIR})
if(NOT TARGET VORO::VORO)
add_library(VORO::VORO UNKNOWN IMPORTED)
set_target_properties(VORO::VORO PROPERTIES
if(NOT TARGET VORO::voro++)
add_library(VORO::voro++ UNKNOWN IMPORTED)
set_target_properties(VORO::voro++ PROPERTIES
IMPORTED_LOCATION "${VORO_LIBRARY}"
INTERFACE_INCLUDE_DIRECTORIES "${VORO_INCLUDE_DIR}")
endif()

1
cmake/Modules/Packages/KOKKOS.cmake
View File

@ -130,6 +130,7 @@ set(KOKKOS_PKG_SOURCES ${KOKKOS_PKG_SOURCES_DIR}/kokkos.cpp
${KOKKOS_PKG_SOURCES_DIR}/atom_vec_kokkos.cpp
${KOKKOS_PKG_SOURCES_DIR}/comm_kokkos.cpp
${KOKKOS_PKG_SOURCES_DIR}/comm_tiled_kokkos.cpp
${KOKKOS_PKG_SOURCES_DIR}/group_kokkos.cpp
${KOKKOS_PKG_SOURCES_DIR}/min_kokkos.cpp
${KOKKOS_PKG_SOURCES_DIR}/min_linesearch_kokkos.cpp
${KOKKOS_PKG_SOURCES_DIR}/neighbor_kokkos.cpp

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

@ -54,5 +54,5 @@ else()
if(NOT VORO_FOUND)
message(FATAL_ERROR "Voro++ library not found. Help CMake to find it by setting VORO_LIBRARY and VORO_INCLUDE_DIR, or set DOWNLOAD_VORO=ON to download it")
endif()
target_link_libraries(lammps PRIVATE VORO::VORO)
target_link_libraries(lammps PRIVATE VORO::voro++)
endif()

9
doc/documentation_conventions.md
View File

@ -10,7 +10,7 @@ Last change: 2022-12-30
In fall 2019, the LAMMPS documentation file format has changed from a
home grown markup designed to generate HTML format files only, to
[reStructuredText](https://docutils.sourceforge.io/rst.html>. For a
[reStructuredText](https://docutils.sourceforge.io/rst.html>). For a
transition period all files in the old .txt format were transparently
converted to .rst and then processed. The `txt2rst tool` is still
included in the distribution to obtain an initial .rst file for legacy
@ -45,8 +45,7 @@ what kind of information and sections are needed.
## Formatting conventions
For headlines we try to follow the conventions posted here:
https://documentation-style-guide-sphinx.readthedocs.io/en/latest/style-guide.html#headings
For headlines we try to follow the conventions posted [here](https://documentation-style-guide-sphinx.readthedocs.io/en/latest/style-guide.html#headings).
It seems to be sufficient to have this consistent only within
any single file and it is not (yet) enforced strictly, but making
this globally consistent makes it easier to move sections around.
@ -64,7 +63,7 @@ Groups of shell commands or LAMMPS input script or C/C++/Python source
code should be typeset into a `.. code-block::` section. A syntax
highlighting extension for LAMMPS input scripts is provided, so `LAMMPS`
can be used to indicate the language in the code block in addition to
`bash`, `c`, `c++`, `console`, `csh`, `diff', `fortran`, `json`, `make`,
`bash`, `c`, `c++`, `console`, `csh`, `diff`, `fortran`, `json`, `make`,
`perl`, `powershell`, `python`, `sh`, or `tcl`, `text`, or `yaml`. When
no syntax style is indicated, no syntax highlighting is performed. When
typesetting commands executed on the shell, please do not prefix
@ -84,7 +83,7 @@ block can be used, followed by multiple `.. tab::` blocks, one
for each alternative. This is only used for HTML output. For other
outputs, the `.. tabs::` directive is transparently removed and
the individual `.. tab::` blocks will be replaced with an
`.. admonition::`` block. Thus in PDF and ePUB output those will
`.. admonition::` block. Thus in PDF and ePUB output those will
be realized as sequential and plain notes.
Special remarks can be highlighted with a `.. note::` block and

2
doc/doxygen/Doxyfile.in
View File

@ -2,7 +2,7 @@
DOXYFILE_ENCODING = UTF-8
PROJECT_NAME = "LAMMPS Programmer's Guide"
PROJECT_NUMBER = "4 May 2022"
PROJECT_NUMBER = "19 November 2024"
PROJECT_BRIEF = "Documentation of the LAMMPS library interface and Python wrapper"
PROJECT_LOGO = lammps-logo.png
CREATE_SUBDIRS = NO

15
doc/github-development-workflow.md
View File

@ -6,7 +6,9 @@ choices the LAMMPS developers have agreed on. Git and GitHub provide the
tools, but do not set policies, so it is up to the developers to come to
an agreement as to how to define and interpret policies. This document
is likely to change as our experiences and needs change, and we try to
adapt it accordingly. Last change 2023-02-10.
adapt it accordingly.
Last change: 2023-02-10
## Table of Contents
@ -72,7 +74,7 @@ be assigned to signal urgency to merge this pull request quickly.
People can be assigned to review a pull request in two ways:
* They can be assigned manually to review a pull request
by the submitter or a LAMMPS developer
by the submitter or a LAMMPS developer.
* They can be automatically assigned, because a developer's GitHub
handle matches a file pattern in the `.github/CODEOWNERS` file,
which associates developers with the code they contributed and
@ -86,9 +88,9 @@ required before merging, in addition to passing all automated
compilation and unit tests. Merging counts as implicit approval, so
does submission of a pull request (by a LAMMPS developer). So the person
doing the merge may not also submit an approving review. The GitHub
feature, that reviews from code owners are "hard" reviews (i.e. they
must all approve before merging is allowed), is currently disabled.
It is in the discretion of the merge maintainer to assess when a
feature that reviews from code owners are "hard" reviews (i.e. they
must all approve before merging is allowed) is currently disabled.
It is at the discretion of the merge maintainer to assess when a
sufficient degree of approval has been reached, especially from external
collaborators. Reviews may be (automatically) dismissed, when the
reviewed code has been changed. Review may be requested a second time.
@ -147,7 +149,8 @@ only contain bug fixes, feature additions to peripheral functionality,
and documentation updates. In between stable releases, bug fixes and
infrastructure updates are back-ported from the "develop" branch to the
"maintenance" branch and occasionally merged into "stable" and published
as update releases.
as update releases. Further explanation of LAMMPS versions can be found
[in the documentation](https://docs.lammps.org/Manual_version.html).
## Project Management

4
doc/lammps.1
View File

@ -1,7 +1,7 @@
.TH LAMMPS "1" "29 August 2024" "2024-08-29"
.TH LAMMPS "1" "19 November 2024" "2024-11-19"
.SH NAME
.B LAMMPS
\- Molecular Dynamics Simulator. Version 29 August 2024
\- Molecular Dynamics Simulator. Version 19 November 2024
.SH SYNOPSIS
.B lmp

9
doc/src/Commands_fix.rst
View File

@ -43,7 +43,7 @@ OPT.
* :doc:`brownian/asphere <fix_brownian>`
* :doc:`brownian/sphere <fix_brownian>`
* :doc:`charge/regulation <fix_charge_regulation>`
* :doc:`cmap <fix_cmap>`
* :doc:`cmap (k) <fix_cmap>`
* :doc:`colvars <fix_colvars>`
* :doc:`controller <fix_controller>`
* :doc:`damping/cundall <fix_damping_cundall>`
@ -134,7 +134,7 @@ OPT.
* :doc:`nve/dot <fix_nve_dot>`
* :doc:`nve/dotc/langevin <fix_nve_dotc_langevin>`
* :doc:`nve/eff <fix_nve_eff>`
* :doc:`nve/limit <fix_nve_limit>`
* :doc:`nve/limit (k) <fix_nve_limit>`
* :doc:`nve/line <fix_nve_line>`
* :doc:`nve/manifold/rattle <fix_nve_manifold_rattle>`
* :doc:`nve/noforce <fix_nve_noforce>`
@ -187,10 +187,11 @@ OPT.
* :doc:`qeq/slater <fix_qeq>`
* :doc:`qmmm <fix_qmmm>`
* :doc:`qtb <fix_qtb>`
* :doc:`qtpie/reaxff <fix_qtpie_reaxff>`
* :doc:`rattle <fix_shake>`
* :doc:`reaxff/bonds (k) <fix_reaxff_bonds>`
* :doc:`reaxff/species (k) <fix_reaxff_species>`
* :doc:`recenter <fix_recenter>`
* :doc:`recenter (k) <fix_recenter>`
* :doc:`restrain <fix_restrain>`
* :doc:`rheo <fix_rheo>`
* :doc:`rheo/oxidation <fix_rheo_oxidation>`
@ -268,7 +269,7 @@ OPT.
* :doc:`wall/piston <fix_wall_piston>`
* :doc:`wall/reflect (k) <fix_wall_reflect>`
* :doc:`wall/reflect/stochastic <fix_wall_reflect_stochastic>`
* :doc:`wall/region <fix_wall_region>`
* :doc:`wall/region (k) <fix_wall_region>`
* :doc:`wall/region/ees <fix_wall_ees>`
* :doc:`wall/srd <fix_wall_srd>`
* :doc:`wall/table <fix_wall>`

16
doc/src/Developer_comm_ops.rst
View File

@ -79,19 +79,19 @@ containing ``double`` values. To correctly store integers that may be
64-bit (bigint, tagint, imageint) in the buffer, you need to use the
:ref:`ubuf union <communication_buffer_coding_with_ubuf>` construct.
The *Fix*, *Compute*, and *Dump* classes can also invoke the same kind
of forward and reverse communication operations using the same *Comm*
class methods. Likewise, the same pack/unpack methods and
The *Fix*, *Bond*, *Compute*, and *Dump* classes can also invoke the
same kind of forward and reverse communication operations using the
same *Comm* class methods. Likewise, the same pack/unpack methods and
comm_forward/comm_reverse variables must be defined by the calling
*Fix*, *Compute*, or *Dump* class.
*Fix*, *Bond*, *Compute*, or *Dump* class.
For *Fix* classes, there is an optional second argument to the
For all of these classes, there is an optional second argument to the
*forward_comm()* and *reverse_comm()* call which can be used when the
fix performs multiple modes of communication, with different numbers
of values per atom. The fix should set the *comm_forward* and
class performs multiple modes of communication, with different numbers
of values per atom. The class should set the *comm_forward* and
*comm_reverse* variables to the maximum value, but can invoke the
communication for a particular mode with a smaller value. For this
to work, the *pack_forward_comm()*, etc methods typically use a class
to work, the *pack_forward_comm()*, etc. methods typically use a class
member variable to choose which values to pack/unpack into/from the
buffer.

23
doc/src/Errors_details.rst
View File

@ -54,3 +54,26 @@ header of a data file (e.g. the number of atoms) is larger than the
number of lines provided (e.g. in the corresponding Atoms section)
and then LAMMPS will continue reading into the next section and that
would have a completely different format.
.. _err0003:
Illegal variable command: expected X arguments but found Y
----------------------------------------------------------
This error indicates that there are the wrong number of arguments for a
specific variable command, but a common reason for that is a variable
expression that has whitespace but is not enclosed in single or double
quotes.
To explain, the LAMMPS input parser reads and processes lines. The
resulting line is broken down into "words". Those are usually
individual commands, labels, names, values separated by whitespace (a
space or tab character). For "words" that may contain whitespace, they
have to be enclosed in single (') or double (") quotes. The parser will
then remove the outermost pair of quotes and then pass that string as
"word" to the variable command.
Thus missing quotes or accidental extra whitespace will lead to the
error shown in the header because the unquoted whitespace will result
in the text being broken into more "words", i.e. the variable expression
being split.

41
doc/src/Howto_bioFF.rst
View File

@ -1,5 +1,5 @@
CHARMM, AMBER, COMPASS, and DREIDING force fields
=================================================
CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields
=======================================================
A compact summary of the concepts, definitions, and properties of
force fields with explicit bonded interactions (like the ones discussed
@ -236,6 +236,40 @@ documentation for the formula it computes.
* :doc:`special_bonds <special_bonds>` dreiding
OPLS
----
OPLS (Optimized Potentials for Liquid Simulations) is a general force
field for atomistic simulation of organic molecules in solvent. It was
developed by the `Jorgensen group
<https://traken.chem.yale.edu/oplsaam.html>`_ at Purdue University and
later at Yale University. Multiple versions of the OPLS parameters
exist for united atom representations (OPLS-UA) and for all-atom
representations (OPLS-AA).
This force field is based on atom types mapped to specific functional
groups in organic and biological molecules. Each atom includes a
static, partial atomic charge reflecting the oxidation state of the
element derived from its bonded neighbors :ref:`(Jorgensen)
<howto-jorgensen>` and computed based on increments determined by the
atom type of the atoms bond to it.
The interaction styles listed below compute force field formulas that
are fully or in part consistent with the OPLS style force fields. See
each command's documentation for the formula it computes. Some are only
compatible with a subset of OPLS interactions.
* :doc:`bond_style <bond_harmonic>` harmonic
* :doc:`angle_style <angle_harmonic>` harmonic
* :doc:`dihedral_style <dihedral_opls>` opls
* :doc:`improper_style <improper_cvff>` cvff
* :doc:`improper_style <improper_fourier>` fourier
* :doc:`improper_style <improper_harmonic>` harmonic
* :doc:`pair_style <pair_lj_cut_coul>` lj/cut/coul/cut
* :doc:`pair_style <pair_lj_cut_coul>` lj/cut/coul/long
* :doc:`pair_modify <pair_modify>` geometric
* :doc:`special_bonds <special_bonds>` lj/coul 0.0 0.0 0.5
----------
.. _Typelabel2:
@ -266,3 +300,6 @@ documentation for the formula it computes.
**(Mayo)** Mayo, Olfason, Goddard III (1990). J Phys Chem, 94, 8897-8909. https://doi.org/10.1021/j100389a010
.. _howto-Jorgensen:
**(Jorgensen)** Jorgensen, Tirado-Rives (1988). J Am Chem Soc, 110, 1657-1666. https://doi.org/10.1021/ja00214a001

14
doc/src/Howto_bpm.rst
View File

@ -5,7 +5,11 @@ The BPM package implements bonded particle models which can be used to
simulate mesoscale solids. Solids are constructed as a collection of
particles, which each represent a coarse-grained region of space much
larger than the atomistic scale. Particles within a solid region are
then connected by a network of bonds to provide solid elasticity.
then connected by a network of bonds to model solid elasticity.
There are many names for methods that are based on similar (or
equivalent) capabilities to those in this package, including, but not
limited to, cohesive beam models, bonded DEMs, lattice spring models,
mass spring models, and lattice particle methods.
Unlike traditional bonds in molecular dynamics, the equilibrium bond
length can vary between bonds. Bonds store the reference state. This
@ -42,7 +46,8 @@ Currently, there are two types of bonds included in the BPM package. The
first bond style, :doc:`bond bpm/spring <bond_bpm_spring>`, only applies
pairwise, central body forces. Point particles must have :doc:`bond atom
style <atom_style>` and may be thought of as nodes in a spring
network. Alternatively, the second bond style, :doc:`bond bpm/rotational
network. An optional multibody term can be used to adjust the network's
Poisson's ratio. Alternatively, the second bond style, :doc:`bond bpm/rotational
<bond_bpm_rotational>`, resolves tangential forces and torques arising
with the shearing, bending, and twisting of the bond due to rotation or
displacement of particles. Particles are similar to those used in the
@ -55,8 +60,9 @@ orientation similar to :doc:`fix nve/asphere <fix_nve_asphere>`.
In addition to bond styles, a new pair style :doc:`pair bpm/spring
<pair_bpm_spring>` was added to accompany the bpm/spring bond
style. This pair style is simply a hookean repulsion with similar
velocity damping as its sister bond style.
style. By default, this pair style is simply a hookean repulsion with
similar velocity damping as its sister bond style, but optional
arguments can be used to modify the force.
----------

67
doc/src/Howto_chunk.rst
View File

@ -58,28 +58,30 @@ chunk ID for an individual atom can also be static (e.g. a molecule
ID), or dynamic (e.g. what spatial bin an atom is in as it moves).
Note that this compute allows the per-atom output of other
:doc:`computes <compute>`, :doc:`fixes <fix>`, and
:doc:`variables <variable>` to be used to define chunk IDs for each
atom. This means you can write your own compute or fix to output a
per-atom quantity to use as chunk ID. See the :doc:`Modify <Modify>`
doc pages for info on how to do this. You can also define a :doc:`per-atom variable <variable>` in the input script that uses a formula to
generate a chunk ID for each atom.
:doc:`computes <compute>`, :doc:`fixes <fix>`, and :doc:`variables
<variable>` to be used to define chunk IDs for each atom. This means
you can write your own compute or fix to output a per-atom quantity to
use as chunk ID. See the :doc:`Modify <Modify>` doc pages for info on
how to do this. You can also define a :doc:`per-atom variable
<variable>` in the input script that uses a formula to generate a chunk
ID for each atom.
Fix ave/chunk command:
----------------------
This fix takes the ID of a :doc:`compute chunk/atom <compute_chunk_atom>` command as input. For each chunk,
it then sums one or more specified per-atom values over the atoms in
each chunk. The per-atom values can be any atom property, such as
velocity, force, charge, potential energy, kinetic energy, stress,
etc. Additional keywords are defined for per-chunk properties like
density and temperature. More generally any per-atom value generated
by other :doc:`computes <compute>`, :doc:`fixes <fix>`, and :doc:`per-atom variables <variable>`, can be summed over atoms in each chunk.
This fix takes the ID of a :doc:`compute chunk/atom
<compute_chunk_atom>` command as input. For each chunk, it then sums
one or more specified per-atom values over the atoms in each chunk. The
per-atom values can be any atom property, such as velocity, force,
charge, potential energy, kinetic energy, stress, etc. Additional
keywords are defined for per-chunk properties like density and
temperature. More generally any per-atom value generated by other
:doc:`computes <compute>`, :doc:`fixes <fix>`, and :doc:`per-atom
variables <variable>`, can be summed over atoms in each chunk.
Similar to other averaging fixes, this fix allows the summed per-chunk
values to be time-averaged in various ways, and output to a file. The
fix produces a global array as output with one row of values per
chunk.
fix produces a global array as output with one row of values per chunk.
Compute \*/chunk commands:
--------------------------
@ -97,17 +99,20 @@ category:
* :doc:`compute torque/chunk <compute_vcm_chunk>`
* :doc:`compute vcm/chunk <compute_vcm_chunk>`
They each take the ID of a :doc:`compute chunk/atom <compute_chunk_atom>` command as input. As their names
indicate, they calculate the center-of-mass, radius of gyration,
moments of inertia, mean-squared displacement, temperature, torque,
and velocity of center-of-mass for each chunk of atoms. The :doc:`compute property/chunk <compute_property_chunk>` command can tally the
count of atoms in each chunk and extract other per-chunk properties.
They each take the ID of a :doc:`compute chunk/atom
<compute_chunk_atom>` command as input. As their names indicate, they
calculate the center-of-mass, radius of gyration, moments of inertia,
mean-squared displacement, temperature, torque, and velocity of
center-of-mass for each chunk of atoms. The :doc:`compute
property/chunk <compute_property_chunk>` command can tally the count of
atoms in each chunk and extract other per-chunk properties.
The reason these various calculations are not part of the :doc:`fix ave/chunk command <fix_ave_chunk>`, is that each requires a more
The reason these various calculations are not part of the :doc:`fix
ave/chunk command <fix_ave_chunk>`, is that each requires a more
complicated operation than simply summing and averaging over per-atom
values in each chunk. For example, many of them require calculation
of a center of mass, which requires summing mass\*position over the
atoms and then dividing by summed mass.
values in each chunk. For example, many of them require calculation of
a center of mass, which requires summing mass\*position over the atoms
and then dividing by summed mass.
All of these computes produce a global vector or global array as
output, with one or more values per chunk. The output can be used in
@ -118,9 +123,10 @@ various ways:
* As input to the :doc:`fix ave/histo <fix_ave_histo>` command to
histogram values across chunks. E.g. a histogram of cluster sizes or
molecule diffusion rates.
* As input to special functions of :doc:`equal-style variables <variable>`, like sum() and max() and ave(). E.g. to
find the largest cluster or fastest diffusing molecule or average
radius-of-gyration of a set of molecules (chunks).
* As input to special functions of :doc:`equal-style variables
<variable>`, like sum() and max() and ave(). E.g. to find the largest
cluster or fastest diffusing molecule or average radius-of-gyration of
a set of molecules (chunks).
Other chunk commands:
---------------------
@ -138,9 +144,10 @@ spatially average per-chunk values calculated by a per-chunk compute.
The :doc:`compute reduce/chunk <compute_reduce_chunk>` command reduces a
peratom value across the atoms in each chunk to produce a value per
chunk. When used with the :doc:`compute chunk/spread/atom <compute_chunk_spread_atom>` command it can
create peratom values that induce a new set of chunks with a second
:doc:`compute chunk/atom <compute_chunk_atom>` command.
chunk. When used with the :doc:`compute chunk/spread/atom
<compute_chunk_spread_atom>` command it can create peratom values that
induce a new set of chunks with a second :doc:`compute chunk/atom
<compute_chunk_atom>` command.
Example calculations with chunks
--------------------------------

17
doc/src/Howto_rheo.rst
View File

@ -15,8 +15,9 @@ details of the system, or develop new capabilities. For instance, the numerics
associated with calculating gradients, reproducing kernels, etc. are separated
into distinct classes to simplify the development of new integration schemes
which can call these calculations. Additional numerical details can be found in
:ref:`(Clemmer) <howto_rheo_clemmer>`. Example movies illustrating some of these
capabilities are found at https://www.lammps.org/movies.html#rheopackage.
:ref:`(Palermo) <howto_rheo_palermo>` and :ref:`(Clemmer) <howto_rheo_clemmer>`.
Example movies illustrating some of these capabilities are found at
https://www.lammps.org/movies.html#rheopackage.
Note, if you simply want to run a traditional SPH simulation, the :ref:`SPH package
<PKG-SPH>` package is likely better suited for your application. It has fewer advanced
@ -70,7 +71,7 @@ particles to solid (e.g. with the :doc:`set <set>` command), (b) create bpm
bonds between the particles (see the :doc:`bpm howto <Howto_bpm>` page for
more details), and (c) use :doc:`pair rheo/solid <pair_rheo_solid>` to
apply repulsive contact forces between distinct solid bodies. Akin to pair rheo,
pair rheo/solid considers a particles fluid/solid phase to determine whether to
pair rheo/solid considers a particle's fluid/solid phase to determine whether to
apply forces. However, unlike pair rheo, pair rheo/solid does obey special bond
settings such that contact forces do not have to be calculated between two bonded
solid particles in the same elastic body.
@ -79,10 +80,10 @@ In systems with thermal evolution, fix rheo/thermal can optionally set a
melting/solidification temperature allowing particles to dynamically swap their
state between fluid and solid when the temperature exceeds or drops below the
critical temperature, respectively. Using the *react* option, one can specify a maximum
bond length and a bond type. Then, when solidifying, particles will search their
bond length and a bond type. Then, when solidifying, particles search their
local neighbors and automatically create bonds with any neighboring solid particles
in range. For BPM bond styles, bonds will then use the immediate position of the two
particles to calculate a reference state. When melting, particles will delete any
in range. For BPM bond styles, bonds then use the immediate position of the two
particles to calculate a reference state. When melting, particles delete any
bonds of the specified type when reverting to a fluid state. Special bonds are updated
as bonds are created/broken.
@ -107,6 +108,10 @@ criteria for creating/deleting a bond or altering force calculations).
----------
.. _howto_rheo_palermo:
**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024).
.. _howto_rheo_clemmer:
**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

4
doc/src/Intro_authors.rst
View File

@ -8,6 +8,8 @@ send an email to all of them at this address: "developers at
lammps.org". General questions about LAMMPS should be posted in the
`LAMMPS forum on MatSci <https://matsci.org/lammps/>`_.
.. We need to keep this file in sync with https://www.lammps.org/authors.html
.. raw:: latex
\small
@ -27,7 +29,7 @@ lammps.org". General questions about LAMMPS should be posted in the
* - `Steve Plimpton <sjp_>`_
- SNL (retired)
- sjplimp at gmail.com
- MD kernels, parallel algorithms & scalability, code structure and design
- original author, MD kernels, parallel algorithms & scalability, code structure and design
* - `Aidan Thompson <at_>`_
- SNL
- athomps at sandia.gov

28
doc/src/Packages_details.rst
View File

@ -2789,14 +2789,15 @@ implements smoothed particle hydrodynamics (SPH) for liquids. See the
related :ref:`MACHDYN package <PKG-MACHDYN>` package for smooth Mach dynamics
(SMD) for solids.
This package contains ideal gas, Lennard-Jones equation of states,
Tait, and full support for complete (i.e. internal-energy dependent)
equations of state. It allows for plain or Monaghans XSPH integration
of the equations of motion. It has options for density continuity or
density summation to propagate the density field. It has
:doc:`set <set>` command options to set the internal energy and density
of particles from the input script and allows the same quantities to
be output with thermodynamic output or to dump files via the :doc:`compute property/atom <compute_property_atom>` command.
This package contains ideal gas, Lennard-Jones equation of states, Tait,
and full support for complete (i.e. internal-energy dependent) equations
of state. It allows for plain or Monaghans XSPH integration of the
equations of motion. It has options for density continuity or density
summation to propagate the density field. It has :doc:`set <set>`
command options to set the internal energy and density of particles from
the input script and allows the same quantities to be output with
thermodynamic output or to dump files via the :doc:`compute
property/atom <compute_property_atom>` command.
**Author:** Georg Ganzenmuller (Fraunhofer-Institute for High-Speed
Dynamics, Ernst Mach Institute, Germany).
@ -2809,6 +2810,17 @@ Dynamics, Ernst Mach Institute, Germany).
* ``examples/PACKAGES/sph``
* https://www.lammps.org/movies.html#sph
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
.. note::
Please also note, that the :ref:`RHEO package <PKG-RHEO>` offers
similar functionality in a more modern and flexible implementation.
----------
.. _PKG-SPIN:

9
doc/src/bond_bpm_rotational.rst
View File

@ -132,8 +132,8 @@ or :doc:`read_restart <read_restart>` commands:
* :math:`k_b` (force*distance/radians units)
* :math:`f_{r,c}` (force units)
* :math:`f_{s,c}` (force units)
* :math:`\tau_{b,c}` (force*distance units)
* :math:`\tau_{t,c}` (force*distance units)
* :math:`\tau_{b,c}` (force*distance units)
* :math:`\gamma_n` (force/velocity units)
* :math:`\gamma_s` (force/velocity units)
* :math:`\gamma_r` (force*distance/velocity units)
@ -154,8 +154,11 @@ page on BPMs.
If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break
during a simulation run. This will prevent some unnecessary calculation.
However, if a bond reaches a damage criterion greater than one,
it will trigger an error.
The recommended bond communication distance no longer depends on bond failure
coefficients (which are ignored) but instead corresponds to the typical heuristic
maximum strain used by typical non-bpm bond styles. Similar behavior to *break no*
can also be attained by setting arbitrarily high values for all four failure
coefficients. One cannot use *break no* with *smooth yes*.
If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed

67
doc/src/bond_bpm_spring.rst
View File

@ -10,7 +10,7 @@ Syntax
bond_style bpm/spring keyword value attribute1 attribute2 ...
* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break*
* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break* or *volume/factor*
.. parsed-literal::
@ -36,6 +36,9 @@ Syntax
*break* value = *yes* or *no*
indicates whether bonds break during a run
*volume/factor* value = *yes* or *no*
indicates whether forces include the volumetric contribution
Examples
""""""""
@ -44,6 +47,9 @@ Examples
bond_style bpm/spring
bond_coeff 1 1.0 0.05 0.1
bond_style bpm/spring volume/factor yes
bond_coeff 1 1.0 0.05 0.1 0.5
bond_style bpm/spring myfix 1000 time id1 id2
dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3]
dump_modify 1 write_header no
@ -97,15 +103,6 @@ approach the critical strain
w = 1.0 - \left( \frac{r - r_0}{r_0 \epsilon_c} \right)^8 .
The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data
<read_data>` or :doc:`read_restart <read_restart>` commands:
* :math:`k` (force/distance units)
* :math:`\epsilon_c` (unit less)
* :math:`\gamma` (force/velocity units)
If the *normalize* keyword is set to *yes*, the elastic bond force will be
normalized by :math:`r_0` such that :math:`k` must be given in force units.
@ -120,8 +117,48 @@ page on BPMs.
If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break
during a simulation run. This will prevent some unnecessary calculation.
However, if a bond reaches a strain greater than :math:`\epsilon_c`,
it will trigger an error.
The recommended bond communication distance no longer depends on the value of
:math:`\epsilon_c` (which is ignored) but instead corresponds to the typical
heuristic maximum strain used by typical non-bpm bond styles. Similar behavior
to *break no* can also be attained by setting an arbitrarily high value of
:math:`\epsilon_c`. One cannot use *break no* with *smooth yes*.
.. versionadded:: TBD
The *volume/factor* keyword toggles whether an additional multibody
contribution is added to he force using the formulation in
:ref:`(Clemmer2) <multibody-Clemmer>`,
.. math::
\alpha_v \left(\left[\frac{V_i + V_j}{V_{0,i} + V_{0,j}}\right]^{1/3} - \frac{r_{ij}}{r_{0,ij}}\right)
where :math:`\alpha_v` is a user specified coefficient and :math:`V_i`
and :math:`V_{0,i}` are estimates of the current and local volume
of atom :math:`i`. These volumes are calculated as the sum of current
or initial bond lengths cubed. In 2D, the volume is replaced with an area
calculated using bond lengths squared and the cube root in the above equation
is accordingly replaced with a square root. This approximation assumes bonds
are evenly distributed on a spherical surface and neglects constant prefactors
which are irrelevant since only the ratio of volumes matters. This term may be
used to adjust the Poisson's ratio.
If a bond is broken (or created), :math:`V_{0,i}` is updated by subtracting
(or adding) that bond's contribution.
The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data
<read_data>` or :doc:`read_restart <read_restart>` commands:
* :math:`k` (force/distance units)
* :math:`\epsilon_c` (unit less)
* :math:`\gamma` (force/velocity units)
Additionally, if *volume/factor* is set to *yes*, a fourth coefficient
must be provided:
* :math:`a_v` (force units)
If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed
@ -213,7 +250,7 @@ Related commands
Default
"""""""
The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes*
The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, *break* = *yes*, and *volume/factor* = *no*
----------
@ -224,3 +261,7 @@ The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *
.. _Groot4:
**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).
.. _multibody-Clemmer:
**(Clemmer2)** Clemmer, Monti, Lechman, Soft Matter, 20, 1702 (2024).

2
doc/src/compute_angle_local.rst
View File

@ -78,7 +78,7 @@ system and output the statistics in various ways:
compute 2 all angle/local eng theta v_cos v_cossq set theta t
dump 1 all local 100 tmp.dump c_1[*] c_2[*]
compute 3 all reduce ave c_2[*]
compute 3 all reduce ave c_2[*] inputs local
thermo_style custom step temp press c_3[*]
fix 10 all ave/histo 10 10 100 -1 1 20 c_2[3] mode vector file tmp.histo

2
doc/src/compute_bond_local.rst
View File

@ -139,7 +139,7 @@ output the statistics in various ways:
compute 2 all bond/local engpot dist v_dsq set dist d
dump 1 all local 100 tmp.dump c_1[*] c_2[*]
compute 3 all reduce ave c_2[*]
compute 3 all reduce ave c_2[*] inputs local
thermo_style custom step temp press c_3[*]
fix 10 all ave/histo 10 10 100 0 6 20 c_2[3] mode vector file tmp.histo

2
doc/src/compute_dihedral_local.rst
View File

@ -76,7 +76,7 @@ angle in the system and output the statistics in various ways:
compute 2 all dihedral/local phi v_cos v_cossq set phi p
dump 1 all local 100 tmp.dump c_1[*] c_2[*]
compute 3 all reduce ave c_2[*]
compute 3 all reduce ave c_2[*] inputs local
thermo_style custom step temp press c_3[*]
fix 10 all ave/histo 10 10 100 -1 1 20 c_2[2] mode vector file tmp.histo

10
doc/src/compute_reduce.rst
View File

@ -206,11 +206,13 @@ IDs and the bond stretch will be printed with thermodynamic output.
The *inputs* keyword allows selection of whether all the inputs are
per-atom or local quantities. As noted above, all the inputs must be
the same kind (per-atom or local). Per-atom is the default setting.
If a compute or fix is specified as an input, it must produce per-atom
or local data to match this setting. If it produces both, e.g. for
the same kind (per-atom or local). Per-atom is the default setting. If
a compute or fix is specified as an input, it must produce per-atom or
local data to match this setting. If it produces both, like for example
the :doc:`compute voronoi/atom <compute_voronoi_atom>` command, then
this keyword selects between them.
this keyword selects between them. If a compute *only* produces local
data, like for example the :doc:`compute bond/local command
<compute_bond_local>`, the setting "inputs local" is *required*.
----------

93
doc/src/compute_reduce_chunk.rst
View File

@ -37,55 +37,57 @@ Description
Define a calculation that reduces one or more per-atom vectors into
per-chunk values. This can be useful for diagnostic output. Or when
used in conjunction with the :doc:`compute chunk/spread/atom <compute_chunk_spread_atom>` command it can be
used to create per-atom values that induce a new set of chunks with a
second :doc:`compute chunk/atom <compute_chunk_atom>` command. An
example is given below.
used in conjunction with the :doc:`compute chunk/spread/atom
<compute_chunk_spread_atom>` command it can be used to create per-atom
values that induce a new set of chunks with a second :doc:`compute
chunk/atom <compute_chunk_atom>` command. An example is given below.
In LAMMPS, chunks are collections of atoms defined by a :doc:`compute chunk/atom <compute_chunk_atom>` command, which assigns each atom
to a single chunk (or no chunk). The ID for this command is specified
as chunkID. For example, a single chunk could be the atoms in a
molecule or atoms in a spatial bin. See the :doc:`compute chunk/atom <compute_chunk_atom>` and :doc:`Howto chunk <Howto_chunk>`
doc pages for details of how chunks can be defined and examples of how
they can be used to measure properties of a system.
In LAMMPS, chunks are collections of atoms defined by a :doc:`compute
chunk/atom <compute_chunk_atom>` command, which assigns each atom to a
single chunk (or no chunk). The ID for this command is specified as
chunkID. For example, a single chunk could be the atoms in a molecule
or atoms in a spatial bin. See the :doc:`compute chunk/atom
<compute_chunk_atom>` and :doc:`Howto chunk <Howto_chunk>` doc pages for
details of how chunks can be defined and examples of how they can be
used to measure properties of a system.
For each atom, this compute accesses its chunk ID from the specified
*chunkID* compute. The per-atom value from an input contributes
to a per-chunk value corresponding the the chunk ID.
*chunkID* compute. The per-atom value from an input contributes to a
per-chunk value corresponding the chunk ID.
The reduction operation is specified by the *mode* setting and is
performed over all the per-atom values from the atoms in each chunk.
The *sum* option adds the pre-atom values to a per-chunk total. The
*min* or *max* options find the minimum or maximum value of the
per-atom values for each chunk.
The *sum* option adds the per-atom values to a per-chunk total. The
*min* or *max* options find the minimum or maximum value of the per-atom
values for each chunk.
Note that only atoms in the specified group contribute to the
reduction operation. If the *chunkID* compute returns a 0 for the
chunk ID of an atom (i.e., the atom is not in a chunk defined by the
:doc:`compute chunk/atom <compute_chunk_atom>` command), that atom will
also not contribute to the reduction operation. An input that is a
compute or fix may define its own group which affects the quantities
it returns. For example, a compute with return a zero value for atoms
that are not in the group specified for that compute.
Note that only atoms in the specified group contribute to the reduction
operation. If the *chunkID* compute returns a 0 for the chunk ID of an
atom (i.e., the atom is not in a chunk defined by the :doc:`compute
chunk/atom <compute_chunk_atom>` command), that atom will also not
contribute to the reduction operation. An input that is a compute or
fix may define its own group which affects the quantities it returns.
For example, a compute will return a zero value for atoms that are not
in the group specified for that compute.
Each listed input is operated on independently. Each input can be the
result of a :doc:`compute <compute>` or :doc:`fix <fix>` or the evaluation
of an atom-style :doc:`variable <variable>`.
result of a :doc:`compute <compute>` or :doc:`fix <fix>` or the
evaluation of an atom-style :doc:`variable <variable>`.
Note that for values from a compute or fix, the bracketed index I can
be specified using a wildcard asterisk with the index to effectively
Note that for values from a compute or fix, the bracketed index I can be
specified using a wildcard asterisk with the index to effectively
specify multiple values. This takes the form "\*" or "\*n" or "m\*" or
"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or the
number of columns in the array (for *mode* = vector), then an asterisk
with no numeric values means all indices from 1 to :math:`N`. A leading
asterisk means all indices from 1 to n (inclusive). A trailing
asterisk means all indices from n to :math:`N` (inclusive). A middle asterisk
means all indices from m to n (inclusive).
"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or
the number of columns in the array (for *mode* = vector), then an
asterisk with no numeric values means all indices from 1 to :math:`N`.
A leading asterisk means all indices from 1 to n (inclusive). A
trailing asterisk means all indices from n to :math:`N` (inclusive). A
middle asterisk means all indices from m to n (inclusive).
Using a wildcard is the same as if the individual columns of the array
had been listed one by one. For example, the following two compute reduce/chunk
commands are equivalent, since the
:doc:`compute property/chunk <compute_property_chunk>` command creates a per-atom
had been listed one by one. For example, the following two compute
reduce/chunk commands are equivalent, since the :doc:`compute
property/chunk <compute_property_chunk>` command creates a per-atom
array with 3 columns:
.. code-block:: LAMMPS
@ -164,13 +166,14 @@ Output info
"""""""""""
This compute calculates a global vector if a single input value is
specified, otherwise a global array is output. The number of columns
in the array is the number of inputs provided. The length of the
vector or the number of vector elements or array rows = the number of
chunks *Nchunk* as calculated by the specified :doc:`compute chunk/atom <compute_chunk_atom>` command. The vector or array can
be accessed by any command that uses global values from a compute as
input. See the :doc:`Howto output <Howto_output>` page for an
overview of LAMMPS output options.
specified, otherwise a global array is output. The number of columns in
the array is the number of inputs provided. The length of the vector or
the number of vector elements or array rows = the number of chunks
*Nchunk* as calculated by the specified :doc:`compute chunk/atom
<compute_chunk_atom>` command. The vector or array can be accessed by
any command that uses global values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
The per-atom values for the vector or each column of the array will be
in whatever :doc:`units <units>` the corresponding input value is in.
@ -183,7 +186,9 @@ Restrictions
Related commands
""""""""""""""""
:doc:`compute chunk/atom <compute_chunk_atom>`, :doc:`compute reduce <compute_reduce>`, :doc:`compute chunk/spread/atom <compute_chunk_spread_atom>`
:doc:`compute chunk/atom <compute_chunk_atom>`,
:doc:`compute reduce <compute_reduce>`,
:doc:`compute chunk/spread/atom <compute_chunk_spread_atom>`
Default
"""""""

6
doc/src/compute_sph_e_atom.rst
View File

@ -33,6 +33,12 @@ particle.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The value of the internal energy will be 0.0 for atoms not in the
specified compute group.

6
doc/src/compute_sph_rho_atom.rst
View File

@ -32,6 +32,12 @@ kernel function interpolation using "pair style sph/rhosum".
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The value of the SPH density will be 0.0 for atoms not in the
specified compute group.

6
doc/src/compute_sph_t_atom.rst
View File

@ -37,6 +37,12 @@ particles, i.e. a Smooth-Particle Hydrodynamics particle.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The value of the internal energy will be 0.0 for atoms not in the
specified compute group.

2
doc/src/compute_stress_curvilinear.rst
View File

@ -87,7 +87,7 @@ This array can be output with :doc:`fix ave/time <fix_ave_time>`,
.. code-block:: LAMMPS
compute 1 all stress/spherical 0 0 0 0.1 10
compute p all stress/spherical 0 0 0 0.1 10
fix 2 all ave/time 100 1 100 c_p[*] file dump_p.out mode vector
The values calculated by this compute are "intensive". The stress

3
doc/src/compute_stress_mop.rst
View File

@ -129,6 +129,9 @@ package <Build_package>` doc page on for more info.
The method is implemented for orthogonal simulation boxes whose
size does not change in time, and axis-aligned planes.
Contributions from bonds, angles, and dihedrals are not compatible
with MPI parallel runs.
The method only works with two-body pair interactions, because it
requires the class method ``Pair::single()`` to be implemented, which is
not possible for manybody potentials. In particular, compute

2
doc/src/compute_temp_partial.rst
View File

@ -46,7 +46,7 @@ degrees of freedom.
A symmetric tensor, stored as a six-element vector, is also calculated
by this compute for use in the computation of a pressure tensor by the
:doc:`compute pressue <compute_pressure>` command. The formula for
:doc:`compute pressure <compute_pressure>` command. The formula for
the components of the tensor is the same as the above expression for
:math:`E_\mathrm{kin}`, except that the 1/2 factor is NOT included and
the :math:`v_i^2` is replaced by :math:`v_{i,x} v_{i,y}` for the

1
doc/src/fix.rst
View File

@ -366,6 +366,7 @@ accelerated styles exist.
* :doc:`qeq/slater <fix_qeq>` - charge equilibration via Slater method
* :doc:`qmmm <fix_qmmm>` - functionality to enable a quantum mechanics/molecular mechanics coupling
* :doc:`qtb <fix_qtb>` - implement quantum thermal bath scheme
* :doc:`qtpie/reaxff <fix_qtpie_reaxff>` - apply QTPIE charge equilibration
* :doc:`rattle <fix_shake>` - RATTLE constraints on bonds and/or angles
* :doc:`reaxff/bonds <fix_reaxff_bonds>` - write out ReaxFF bond information
* :doc:`reaxff/species <fix_reaxff_species>` - write out ReaxFF molecule information

6
doc/src/fix_acks2_reaxff.rst
View File

@ -111,7 +111,8 @@ LAMMPS was built with that package. See the :doc:`Build package
This fix does not correctly handle interactions involving multiple
periodic images of the same atom. Hence, it should not be used for
periodic cell dimensions less than :math:`10~\AA`.
periodic cell dimensions smaller than the non-bonded cutoff radius,
which is typically :math:`10~\AA` for ReaxFF simulations.
This fix may be used in combination with :doc:`fix efield <fix_efield>`
and will apply the external electric field during charge equilibration,
@ -122,7 +123,8 @@ components in non-periodic directions.
Related commands
""""""""""""""""
:doc:`pair_style reaxff <pair_reaxff>`, :doc:`fix qeq/reaxff <fix_qeq_reaxff>`
:doc:`pair_style reaxff <pair_reaxff>`, :doc:`fix qeq/reaxff <fix_qeq_reaxff>`,
:doc:`fix qtpi/reaxff <fix_qtpie_reaxff>`
Default
"""""""

72
doc/src/fix_adapt.rst
View File

@ -119,6 +119,14 @@ style supports it. Note that the :doc:`pair_style <pair_style>` and
to specify these parameters initially; the fix adapt command simply
overrides the parameters.
.. note::
Pair_coeff settings must be made **explicitly** in order for fix
adapt to be able to change them. Settings inferred from mixing
are not suitable. If necessary all mixed settings can be output
to a file using the :doc:`write_coeff command <write_coeff>` and
then the desired mixed pair_coeff settings copied from that file.
The *pstyle* argument is the name of the pair style. If
:doc:`pair_style hybrid or hybrid/overlay <pair_hybrid>` is used,
*pstyle* should be a sub-style name. If there are multiple
@ -323,7 +331,7 @@ 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:`class2 <bond_class2>` | k2,k3,k4,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
| :doc:`fene <bond_fene>` | k,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
@ -331,6 +339,8 @@ sub-style name. The bond styles that currently work with fix adapt are:
+-----------------------------------------------------+---------------------------+------------+
| :doc:`fene/nm <bond_fene>` | k,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
| :doc:`gaussian <bond_gaussian>` | alpha,width,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
| :doc:`gromos <bond_gromos>` | k,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
| :doc:`harmonic <bond_harmonic>` | k,r0 | type bonds |
@ -343,9 +353,9 @@ sub-style name. The bond styles that currently work with fix adapt are:
+-----------------------------------------------------+---------------------------+------------+
| :doc:`mm3 <bond_mm3>` | k,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
| :doc:`morse <bond_morse>` | r0 | type bonds |
| :doc:`morse <bond_morse>` | d0,alpha,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
| :doc:`nonlinear <bond_nonlinear>` | epsilon,r0 | type bonds |
| :doc:`nonlinear <bond_nonlinear>` | lamda,epsilon,r0 | type bonds |
+-----------------------------------------------------+---------------------------+------------+
----------
@ -369,31 +379,37 @@ all types from 1 to :math:`N`. A leading asterisk means all types from
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:`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 |
+--------------------------------------------------------------------+-----------------+-------------+
+--------------------------------------------------------------------+--------------------+-------------+
| :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/delta <angle_cosine_delta>` | k | type angles |
+--------------------------------------------------------------------+--------------------+-------------+
| :doc:`cosine/periodic <angle_cosine_periodic>` | k,b,n | type angles |
+--------------------------------------------------------------------+--------------------+-------------+
| :doc:`cosine/squared <angle_cosine_squared>` | k,theta0 | 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:`gaussian <angle_gaussian>` | alpha,width,theta0 | 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.

22
doc/src/fix_adapt_fep.rst
View File

@ -116,12 +116,22 @@ style supports it. Note that the :doc:`pair_style <pair_style>` and
to specify these parameters initially; the fix adapt command simply
overrides the parameters.
The *pstyle* argument is the name of the pair style. If :doc:`pair_style hybrid or hybrid/overlay <pair_hybrid>` is used, *pstyle* should be
a sub-style name. For example, *pstyle* could be specified as "soft"
or "lubricate". The *pparam* argument is the name of the parameter to
change. This is the current list of pair styles and parameters that
can be varied by this fix. See the doc pages for individual pair
styles and their energy formulas for the meaning of these parameters:
.. note::
Pair_coeff settings must be made **explicitly** in order for fix
adapt/fep to be able to change them. Settings inferred from mixing
are not suitable. If necessary all mixed settings can be output
to a file using the :doc:`write_coeff command <write_coeff>` and
then the desired mixed pair_coeff settings copied from that file.
The *pstyle* argument is the name of the pair style. If
:doc:`pair_style hybrid or hybrid/overlay <pair_hybrid>` is used,
*pstyle* should be a sub-style name. For example, *pstyle* could be
specified as "soft" or "lubricate". The *pparam* argument is the name
of the parameter to change. This is the current list of pair styles and
parameters that can be varied by this fix. See the doc pages for
individual pair styles and their energy formulas for the meaning of
these parameters:
+------------------------------------------------------------------------------+-------------------------+------------+
| :doc:`born <pair_born>` | a,b,c | type pairs |

13
doc/src/fix_atom_swap.rst
View File

@ -119,15 +119,14 @@ groups of atoms that have different charges, these charges will not be
changed when the atom types change.
Since this fix computes total potential energies before and after
proposed swaps, so even complicated potential energy calculations are
OK, including the following:
proposed swaps, even complicated potential energy calculations are
acceptable, including the following:
* long-range electrostatics (:math:`k`-space)
* many body pair styles
* hybrid pair styles
* eam pair styles
* hybrid pair styles (with restrictions)
* EAM pair styles
* triclinic systems
* need to include potential energy contributions from other fixes
Some fixes have an associated potential energy. Examples of such fixes
include: :doc:`efield <fix_efield>`, :doc:`gravity <fix_gravity>`,
@ -181,6 +180,10 @@ This fix is part of the MC package. It is only enabled if LAMMPS was
built with that package. See the :doc:`Build package <Build_package>`
doc page for more info.
When this fix is used with a :doc:`hybrid pair style <pair_hybrid>`
system, only swaps between atom types of the same sub-style (or
combination of sub-styles) are permitted.
This fix cannot be used with systems that do not have per-type masses
(e.g. atom style sphere) since the implemented algorithm pre-computes
velocity rescaling factors from per-type masses and ignores any per-atom

9
doc/src/fix_cmap.rst
View File

@ -1,8 +1,11 @@
.. index:: fix cmap
.. index:: fix cmap/kk
fix cmap command
================
Accelerator Variants: *cmap/kk*
Syntax
""""""
@ -141,6 +144,12 @@ outermost level.
MUST not disable the :doc:`fix_modify <fix_modify>` *energy* option
for this fix.
----------
.. include:: accel_styles.rst
----------
Restrictions
""""""""""""

13
doc/src/fix_colvars.rst
View File

@ -1,8 +1,11 @@
.. index:: fix colvars
.. index:: fix colvars/kk
fix colvars command
===================
Accelerator Variants: *colvars/kk*
Syntax
""""""
@ -118,6 +121,16 @@ thermostat target temperature.
The *seed* keyword contains the seed for the random number generator
that will be used in the colvars module.
----------
.. note::
Fix colvars/kk is not really ported to KOKKOS, since the colvars
library has not been ported to KOKKOS. It merely has some
optimizations to reduce the data transfers between host and device
for KOKKOS with GPUs.
----------
Restarting
""""""""""

22
doc/src/fix_external.rst
View File

@ -139,9 +139,9 @@ properties are:
void set_vector_length(int n);
void set_vector(int idx, double val);
These allow to set per-atom energy contributions, per-atom stress
contributions, the length and individual values of a global vector
of properties that the caller code may want to communicate to LAMMPS
These enable setting per-atom energy and per-atom stress contributions,
the length and individual values of a global vector of properties that
the caller code may want to communicate to LAMMPS
(e.g. for use in :doc:`fix ave/time <fix_ave_time>` or in
:doc:`equal-style variables <variable>` or for
:doc:`custom thermo output <thermo_style>`.
@ -173,9 +173,19 @@ stress/atom <compute_stress_atom>` commands. The former can be
accessed by :doc:`thermodynamic output <thermo_style>`. The default
setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.
This fix computes a global scalar which can be accessed by various
:doc:`output commands <Howto_output>`. The scalar is the potential
energy discussed above. The scalar stored by this fix is "extensive".
This fix computes a global scalar, a global vector, and a per-atom array
which can be accessed by various :doc:`output commands <Howto_output>`.
The scalar is the potential energy discussed above. The scalar stored
by this fix is "extensive".
The global vector has a custom length and needs to be set by the external
program using the
:cpp:func:`lammps_fix_external_set_vector() <lammps_fix_external_set_vector>`
and :cpp:func:`lammps_fix_external_set_vector_length()
<lammps_fix_external_set_vector_length>`
calls of the LAMMPS library interface or the equivalent call of the Python
or Fortran modules.
The per-atom array has 3 values for each atom and is the applied external
force.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command.

2
doc/src/fix_halt.rst
View File

@ -101,7 +101,7 @@ hstyle = bondmax option.
.. code-block:: LAMMPS
compute bdist all bond/local dist
compute bmax all reduce max c_bdist
compute bmax all reduce max c_bdist inputs local
variable bondmax equal c_bmax
Thus these two versions of a fix halt command will do the same thing:

6
doc/src/fix_langevin.rst
View File

@ -231,12 +231,6 @@ the particles. As described below, this energy can then be printed
out or added to the potential energy of the system to monitor energy
conservation.
.. note::
This accumulated energy does NOT include kinetic energy removed
by the *zero* flag. LAMMPS will print a warning when both options are
active.
The keyword *zero* can be used to eliminate drift due to the
thermostat. Because the random forces on different atoms are
independent, they do not sum exactly to zero. As a result, this fix

9
doc/src/fix_nve_limit.rst
View File

@ -1,8 +1,11 @@
.. index:: fix nve/limit
.. index:: fix nve/limit/kk
fix nve/limit command
=====================
Accelerator Variants: *nve/limit/kk*
Syntax
""""""
@ -79,6 +82,12 @@ is "extensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
----------
.. include:: accel_styles.rst
----------
Restrictions
""""""""""""
none

2
doc/src/fix_qeq.rst
View File

@ -190,7 +190,7 @@ on atoms via the matrix inversion method. A tolerance of 1.0e-6 is
usually a good number. Keyword *alpha* can be used to change the Slater
type orbital exponent.
.. versionadded:: TBD
.. versionadded:: 19Nov2024
The *qeq/ctip* style describes partial charges on atoms in the same way
as style *qeq/shielded* but also enables the definition of charge

6
doc/src/fix_qeq_reaxff.rst
View File

@ -124,7 +124,8 @@ LAMMPS was built with that package. See the :doc:`Build package
This fix does not correctly handle interactions involving multiple
periodic images of the same atom. Hence, it should not be used for
periodic cell dimensions less than 10 Angstroms.
periodic cell dimensions smaller than the non-bonded cutoff radius,
which is typically :math:`10~\AA` for ReaxFF simulations.
This fix may be used in combination with :doc:`fix efield <fix_efield>`
and will apply the external electric field during charge equilibration,
@ -138,7 +139,8 @@ as an atom-style variable using the *potential* keyword for `fix efield`.
Related commands
""""""""""""""""
:doc:`pair_style reaxff <pair_reaxff>`, :doc:`fix qeq/shielded <fix_qeq>`
:doc:`pair_style reaxff <pair_reaxff>`, :doc:`fix qeq/shielded <fix_qeq>`,
:doc:`fix acks2/reaxff <fix_acks2_reaxff>`, :doc:`fix qtpie/reaxff <fix_qtpie_reaxff>`
Default
"""""""

200
doc/src/fix_qtpie_reaxff.rst Normal file
View File

@ -0,0 +1,200 @@
.. index:: fix qtpie/reaxff
fix qtpie/reaxff command
========================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID qtpie/reaxff Nevery cutlo cuthi tolerance params gfile args
* ID, group-ID are documented in :doc:`fix <fix>` command
* qtpie/reaxff = style name of this fix command
* Nevery = perform QTPIE every this many steps
* cutlo,cuthi = lo and hi cutoff for Taper radius
* tolerance = precision to which charges will be equilibrated
* params = reaxff or a filename
* gfile = the name of a file containing Gaussian orbital exponents
* one or more keywords or keyword/value pairs may be appended
.. parsed-literal::
keyword = *maxiter*
*maxiter* N = limit the number of iterations to *N*
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff exp.qtpie
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 params.qtpie exp.qtpie maxiter 500
Description
"""""""""""
.. versionadded:: 19Nov2024
The QTPIE charge equilibration method is an extension of the QEq charge
equilibration method. With QTPIE, the partial charges on individual atoms
are computed by minimizing the electrostatic energy of the system in the
same way as the QEq method but where the absolute electronegativity,
:math:`\chi_i`, of each atom in the QEq charge equilibration scheme
:ref:`(Rappe and Goddard) <Rappe3>` is replaced with an effective
electronegativity given by :ref:`(Chen) <qtpie-Chen>`
.. math::
\chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j) S_{ij}}
{\sum_{m=1}^{N}S_{im}},
which acts to penalize long-range charge transfer seen with the QEq charge
equilibration scheme. In this equation, :math:`N` is the number of atoms in
the system and :math:`S_{ij}` is the overlap integral between atom :math:`i`
and atom :math:`j`.
The effect of an external electric field can be incorporated into the QTPIE
method by modifying the absolute or effective electronegativities of each
atom :ref:`(Chen) <qtpie-Chen>`. This fix models the effect of an external
electric field by using the effective electronegativity given in
:ref:`(Gergs) <Gergs>`:
.. math::
\chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_i - \phi_j) S_{ij}}
{\sum_{m=1}^{N}S_{im}},
where :math:`\phi_i` and :math:`\phi_j` are the electric
potentials at the positions of atom :math:`i` and :math:`j`
due to the external electric field.
This fix is typically used in conjunction with the ReaxFF force
field model as implemented in the :doc:`pair_style reaxff <pair_reaxff>`
command, but it can be used with any potential in LAMMPS, so long as it
defines and uses charges on each atom. For more technical details about the
charge equilibration performed by `fix qtpie/reaxff`, which is the same as in
:doc:`fix qeq/reaxff <fix_qeq_reaxff>` except for the use of
:math:`\chi_{\mathrm{eff},i}`, please refer to :ref:`(Aktulga) <qeq-Aktulga2>`.
To be explicit, this fix replaces :math:`\chi_k` of eq. 3 in
:ref:`(Aktulga) <qeq-Aktulga2>` with :math:`\chi_{\mathrm{eff},k}`.
This fix requires the absolute electronegativity, :math:`\chi`, in eV, the
self-Coulomb potential, :math:`\eta`, in eV, and the shielded Coulomb
constant, :math:`\gamma`, in :math:`\AA^{-1}`. If the *params* setting above
is the word "reaxff", then these are extracted from the
:doc:`pair_style reaxff <pair_reaxff>` command and the ReaxFF force field
file it reads in. If a file name is specified for *params*, then the
parameters are taken from the specified file and the file must contain
one line for each atom type. The latter form must be used when performing
QTPIE with a non-ReaxFF potential. Each line should be formatted as follows,
ensuring that the parameters are given in units of eV, eV, and :math:`\AA^{-1}`,
respectively:
.. parsed-literal::
itype chi eta gamma
where *itype* is the atom type from 1 to Ntypes. Note that eta is
defined here as twice the eta value in the ReaxFF file.
The overlap integrals in the equation for :math:`\chi_{\mathrm{eff},i}`
are computed by using normalized 1s Gaussian type orbitals. The Gaussian
orbital exponents, :math:`\alpha`, that are needed to compute the overlap
integrals are taken from the file given by *gfile*.
This file must contain one line for each atom type and provide the Gaussian
orbital exponent for each atom type in units of inverse square Bohr radius.
Each line should be formatted as follows:
.. parsed-literal::
itype alpha
Empty lines or any text following the pound sign (#) are ignored. An example
*gfile* for a system with two atom types is
.. parsed-literal::
# An example gfile. Exponents are taken from Table 2.2 of Chen, J. (2009).
# Theory and applications of fluctuating-charge models.
# The units of the exponents are 1 / (Bohr radius)^2 .
1 0.2240 # O
2 0.5434 # H
The optional *maxiter* keyword allows changing the max number
of iterations in the linear solver. The default value is 200.
.. note::
In order to solve the self-consistent equations for electronegativity
equalization, LAMMPS imposes the additional constraint that all the
charges in the fix group must add up to zero. The initial charge
assignments should also satisfy this constraint. LAMMPS will print a
warning if that is not the case.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files
<restart>`. This fix computes a global scalar (the number of
iterations) and a per-atom vector (the effective electronegativity), which
can be accessed by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command.
This fix is invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix is part of the REAXFF package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This fix does not correctly handle interactions involving multiple
periodic images of the same atom. Hence, it should not be used for
periodic cell dimensions smaller than the non-bonded cutoff radius,
which is typically :math:`10~\AA` for ReaxFF simulations.
This fix may be used in combination with :doc:`fix efield <fix_efield>`
and will apply the external electric field during charge equilibration,
but there may be only one fix efield instance used and the electric field
must be applied to all atoms in the system. Consequently, `fix efield` must
be used with *group-ID* all and must not be used with the keyword *region*.
Equal-style variables can be used for electric field vector
components without any further settings. Atom-style variables can be used
for spatially-varying electric field vector components, but the resulting
electric potential must be specified as an atom-style variable using
the *potential* keyword for `fix efield`.
Related commands
""""""""""""""""
:doc:`pair_style reaxff <pair_reaxff>`, :doc:`fix qeq/reaxff <fix_qeq_reaxff>`,
:doc:`fix acks2/reaxff <fix_acks2_reaxff>`
Default
"""""""
maxiter 200
----------
.. _Rappe3:
**(Rappe)** Rappe and Goddard III, Journal of Physical Chemistry, 95,
3358-3363 (1991).
.. _qtpie-Chen:
**(Chen)** Chen, Jiahao. Theory and applications of fluctuating-charge models.
University of Illinois at Urbana-Champaign, 2009.
.. _Gergs:
**(Gergs)** Gergs, Dirkmann and Mussenbrock.
Journal of Applied Physics 123.24 (2018).
.. _qeq-Aktulga2:
**(Aktulga)** Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38,
245-259 (2012).

9
doc/src/fix_recenter.rst
View File

@ -1,8 +1,11 @@
.. index:: fix recenter
.. index:: fix recenter/kk
fix recenter command
====================
Accelerator Variants: *recenter/kk*
Syntax
""""""
@ -113,6 +116,12 @@ The scalar and vector values calculated by this fix are "extensive".
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
----------
.. include:: accel_styles.rst
----------
Restrictions
""""""""""""

113
doc/src/fix_rheo.rst
View File

@ -16,21 +16,36 @@ Syntax
* kstyle = *quintic* or *RK0* or *RK1* or *RK2*
* zmin = minimal number of neighbors for reproducing kernels
* zero or more keyword/value pairs may be appended to args
* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or *shift* or *rho/sum* or *density* or *self/mass* or *speed/sound*
* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or *shift* or *rho/sum* or *density* or *speed/sound*
.. parsed-literal::
*thermal* values = none, turns on thermal evolution
*interface/reconstruct* values = none, reconstructs interfaces with solid particles
*surface/detection* values = *sdstyle* *limit* *limit/splash*
*sdstyle* = *coordination* or *divergence*
*limit* = threshold for surface particles
*limit/splash* = threshold for splash particles
*shift* values = none, turns on velocity shifting
*rho/sum* values = none, uses the kernel to compute the density of particles
*self/mass* values = none, a particle uses its own mass in a rho summation
*density* values = *rho01*, ... *rho0N* (density)
*speed/sound* values = *cs0*, ... *csN* (velocity)
*thermal* turns on thermal evolution
values = none
*interface/reconstruct* reconstructs interfaces with solid particles
values = none
*surface/detection* detects free-surfaces with an absence of particles
values = *sdstyle* *limit* *limit/splash*
*sdstyle* = *coordination* or *divergence*
*limit* = threshold for surface particles
*limit/splash* = threshold for splash particles (unitless)
*shift* turns on velocity shifting
values = none
optional args = *exclude/type* or *scale/cross/type*
*exclude/type* values = *types*
*types* = list of types
*scale/cross/type* values = *shiftscale* *cmin* *wmin*
*shiftscale* = fraction of shifting in normal direction to preserve (unitless)
*cmin* = minimum color function value required for scaling (unitless)
*wmin* = minimum local same-type support required for any shifting (unitless)
*rho/sum* density evolution performed by a kernel summation
values = none
optional args = *self/mass*
*self/mass* values = none, a particle uses its own mass in summation
*density* specify equilibrium densities for each atom type
values = *rho01*, ... *rho0N* (density)
*speed/sound* specify speeds of sound for each atom type
values = *cs0*, ... *csN* (velocity)
Examples
""""""""
@ -39,6 +54,8 @@ Examples
fix 1 all rheo 3.0 quintic 0 thermal density 0.1 0.1 speed/sound 10.0 1.0
fix 1 all rheo 3.0 RK1 10 shift surface/detection coordination 40
fix 1 all rheo 3.0 RK1 10 shift exclude/type 2*4 scale/cross/type 0.05 0.02 0.5
fix 1 all rheo 3.0 RK1 10 rhosum self/mass
Description
"""""""""""
@ -46,8 +63,10 @@ Description
.. versionadded:: 29Aug2024
Perform time integration for RHEO particles, updating positions, velocities,
and densities. For an overview of other features available in the RHEO package,
see :doc:`the RHEO howto <Howto_rheo>`.
and densities. For a detailed breakdown of the integration timestep and
numerical details, see :ref:`(Palermo) <fix_rheo_palermo>`. For an overview
and list of other features available in the RHEO package, see
:doc:`the RHEO howto <Howto_rheo>`.
The type of kernel is specified using *kstyle* and the cutoff is *cut*. Four
kernels are currently available. The *quintic* kernel is a standard quintic
@ -70,16 +89,51 @@ and velocity of solid particles are alternatively reconstructed for every
fluid-solid interaction to ensure no-slip and pressure-balanced boundaries.
This is done by estimating the location of the fluid-solid interface and
extrapolating fluid particle properties across the interface to calculate a
temporary apparent density and velocity for a solid particle.
temporary apparent density and velocity for a solid particle. The numerical
details are the same as those described in
:ref:`(Palermo) <fix_rheo_palermo>` except there is an additional
restriction that the reconstructed solid density cannot be less than the
equilibrium density. This prevents fluid particles from sticking to solid
surfaces.
A modified form of Fickian particle shifting can be enabled with the
*shift* keyword. This effectively shifts particle positions to generate a
more uniform spatial distribution. Shifting currently does not consider the
more uniform spatial distribution. By default, shifting does not consider the
type of a particle and therefore may be inappropriate in systems consisting
of multiple fluid phases.
of multiple atom types representing multiple fluid phases. However, two
optional sub-arguments can follow the *shift* keyword, *exclude/type* and
*scale/cross/type* to adjust shifting at fluid interfaces.
In systems with free surfaces, the *surface/detection* keyword can be used
to classify the location of particles as being within the bulk fluid, on a
The *exclude/type* option lets the user specify a list of atom types which
are not shifted, *types*. A wild-card asterisk can be used in place
of or in conjunction with the *types* argument to toggle shifting for
multiple atom types. This takes the form "\*" or "\*n" or "m\*"
or "m\*n". If :math:`N` is the number of atom types, then an asterisk with
no numeric values means all types from 1 to :math:`N`. A leading asterisk
means all types from 1 to n (inclusive). A trailing asterisk means all types
from m to :math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
The *scale/cross/type* option is designed to handle interfaces between fluids
made up of different atom types. Similar to the method by
:ref:`(Yang) <fix_rheo_yang>`, a color function is calculated and used to
estimate a local interfacial normal vector. Shifting along this normal direction
is rescaled by a factor of *scaleshift*, such that a value of *scaleshift* of
zero implies there is no shifting in the normal direction and a value of
*scaleshift* of one implies no change in behavior. This scaling is only applied
to atoms with a color function value greater than *cmin*. To handle scenarios
of a small inclusion of one fluid type (e.g. a single atom) inside another,
the degree of same-type support is calculated
.. math::
W_{i,\mathrm{same}} = \sum_{j} W_{ij} \delta_{ij}
where :math:`\delta_{ij}` is zero if atoms :math:`i` and :math:`j` have different
types but unity otherwise. If :math:`W_{i,\mathrm{same}}` is ever less than the
specified value of *wmin*, shifting is turned off for particle :math:`i`
In systems with free surfaces (atom-vacuum), the *surface/detection* keyword
can classify the location of particles as being within the bulk fluid, on a
free surface, or isolated from other particles in a splash or droplet.
Shifting is then disabled in the normal direction away from the free surface
to prevent particles from diffusing away. Surface detection can also be used
@ -101,10 +155,9 @@ threshold for this classification is set by the numerical value of
By default, RHEO integrates particles' densities using a mass diffusion
equation. Alternatively, one can update densities every timestep by performing
a kernel summation of the masses of neighboring particles by specifying the *rho/sum*
keyword.
The *self/mass* keyword modifies the behavior of the density summation in *rho/sum*.
Typically, the density :math:`\rho` of a particle is calculated as the sum over neighbors
keyword. Following this keyword, one may include the optional *self/mass* sub-argument
which modifies the behavior of the density summation. Typically, the density
:math:`\rho` of a particle is calculated as the sum over neighbors
.. math::
\rho_i = \sum_{j} W_{ij} M_j
@ -120,7 +173,9 @@ equilibrium density *rho0*.
The *speed/sound* keyword is used to specify the speed of sound of each of the
N particle types. It must be followed by N numerical values specifying each type's
speed of sound *cs*.
speed of sound *cs*. These values may be ignored if the pressure equation of
state has a non-constant speed of sound, as discussed further in
:doc:`fix rheo/pressure <fix_rheo_pressure>`.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@ -163,6 +218,14 @@ Default
----------
.. _fix_rheo_palermo:
**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, Phys. Fluids, 36, 113337 (2024).
.. _fix_rheo_yang:
**(Yang)** Yang, Rakhsha, Hu, Negrut, J. Comp. Physics, 458, 111079 (2022).
.. _fix_rheo_hu:
**(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006).
**(Hu)** Hu, and Adams, J. Comp. Physics, 213, 844-861 (2006).

52
doc/src/fix_rheo_pressure.rst
View File

@ -14,13 +14,16 @@ Syntax
* rheo/pressure = style name of this fix command
* one or more types and pressure styles must be appended
* types = lists of types (see below)
* pstyle = *linear* or *taitwater* or *cubic*
* pstyle = *linear* or *tait/water* or *tait/general* or *cubic* or *ideal/gas* or *background*
.. parsed-literal::
*linear* args = none
*taitwater* args = none
*tait/water* args = none
*tait/general* args = exponent :math:`gamma` (unitless)
*cubic* args = cubic prefactor :math:`A_3` (pressure/density\^2)
*ideal/gas* args = heat capacity ratio :math:`gamma` (unitless)
*background* args = background pressure :math:`P[b]` (pressure)
Examples
""""""""
@ -29,6 +32,7 @@ Examples
fix 1 all rheo/pressure * linear
fix 1 all rheo/pressure 1 linear 2 cubic 10.0
fix 1 all rheo/pressure * linear * background 0.1
Description
"""""""""""
@ -40,13 +44,12 @@ define different equations of state for different atom types. An equation
must be specified for every atom type.
One first defines the atom *types*. A wild-card asterisk can be used in place
of or in conjunction with the *types* argument to set the coefficients for
multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*"
or "m\*n". If :math:`N` is the number of atom types, then an asterisk with
no numeric values means all types from 1 to :math:`N`. A leading asterisk
means all types from 1 to n (inclusive). A trailing asterisk means all types
from m to :math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
of or in conjunction with the *types* argument to set values for multiple atom
types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is
the number of atom types, then an asterisk with no numeric values means all types
from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive).
A trailing asterisk means all types from m to :math:`N` (inclusive). A middle
asterisk means all types from m to n (inclusive).
The *types* definition is followed by the pressure style, *pstyle*. Current
options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear
@ -54,7 +57,7 @@ equation of state with a particle pressure :math:`P` calculated as
.. math::
P = c (\rho - \rho_0)
P = c^2 (\rho - \rho_0)
where :math:`c` is the speed of sound, :math:`\rho_0` is the equilibrium density,
and :math:`\rho` is the current density of a particle. The numerical values of
@ -63,14 +66,39 @@ is a cubic equation of state which has an extra argument :math:`A_3`,
.. math::
P = c ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) .
P = c^2 ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) .
Style *taitwater* is Tait's equation of state:
Style *tait/water* is Tait's equation of state:
.. math::
P = \frac{c^2 \rho_0}{7} \biggl[\left(\frac{\rho}{\rho_0}\right)^{7} - 1\biggr].
Style *tait/general* generalizes this equation of state
.. math::
P = \frac{c^2 \rho_0}{\gamma} \biggl[\left(\frac{\rho}{\rho_0}\right)^{\gamma} - 1\biggr].
where :math:`\gamma` is an exponent.
Style *ideal/gas* is the ideal gas equation of state
.. math::
P = (\gamma - 1) \rho e
where :math:`\gamma` is the heat capacity ratio and :math:`e` is the internal energy of
a particle per unit mass. This style is only compatible with atom style rheo/thermal.
Note that when using this style, the speed of sound is no longer constant such that the
value of :math:`c` specified in :doc:`fix rheo <fix_rheo>` is not used.
The *background* style acts differently than the rest as it
only adds a constant background pressure shift :math:`P[b]`
to all atoms of the designated types. Therefore, this style
must be used in conjunction with another style that specifies
an equation of state.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

14
doc/src/fix_rheo_thermal.rst
View File

@ -70,13 +70,13 @@ of the energy is used to shift energies. This may be inappropriate in systems
with multiple atom types with different specific heats.
For each property, one must first define a list of atom types. A wild-card
asterisk can be used in place of or in conjunction with the *types* argument
to set the coefficients for multiple pairs of atom types. This takes the
form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is the number of atom
types, then an asterisk with no numeric values means all types from 1 to
:math:`N`. A leading asterisk means all types from 1 to n (inclusive).
A trailing asterisk means all types from m to :math:`N` (inclusive). A
middle asterisk means all types from m to n (inclusive).
asterisk can be used in place of or in conjunction with the *types* argument to
set values for multiple atom types. This takes the form "\*" or "\*n" or "m\*"
or "m\*n". If :math:`N` is the number of atom types, then an asterisk with no
numeric values means all types from 1 to :math:`N`. A leading asterisk means
all types from 1 to n (inclusive). A trailing asterisk means all types from m
to :math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
The *types* definition for each property is followed by the style. Currently,
the only option is *constant*. Style *constant* simply applies a constant value

13
doc/src/fix_rheo_viscosity.rst
View File

@ -45,13 +45,12 @@ viscosities for different atom types, but a viscosity must be specified for
every atom type.
One first defines the atom *types*. A wild-card asterisk can be used in place
of or in conjunction with the *types* argument to set the coefficients for
multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*"
or "m\*n". If :math:`N` is the number of atom types, then an asterisk with
no numeric values means all types from 1 to :math:`N`. A leading asterisk
means all types from 1 to n (inclusive). A trailing asterisk means all types
from m to :math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
of or in conjunction with the *types* argument to set values for multiple atom
types. This takes the form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is
the number of atom types, then an asterisk with no numeric values means all types
from 1 to :math:`N`. A leading asterisk means all types from 1 to n (inclusive).
A trailing asterisk means all types from m to :math:`N` (inclusive). A middle
asterisk means all types from m to n (inclusive).
The *types* definition is followed by the viscosity style, *vstyle*. Two
options are available, *constant* and *power*. Style *constant* simply

6
doc/src/fix_sph.rst
View File

@ -32,6 +32,12 @@ Hydrodynamics.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

6
doc/src/fix_sph_stationary.rst
View File

@ -32,6 +32,12 @@ space. SPH stands for Smoothed Particle Hydrodynamics.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

23
doc/src/fix_spring_self.rst
View File

@ -13,7 +13,7 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* spring/self = style name of this fix command
* K = spring constant (force/distance units)
* K = spring constant (force/distance units), can be a variable (see below)
* dir = xyz, xy, xz, yz, x, y, or z (optional, default: xyz)
Examples
@ -22,6 +22,7 @@ Examples
.. code-block:: LAMMPS
fix tether boundary-atoms spring/self 10.0
fix var all spring/self v_kvar
fix zrest move spring/self 10.0 z
Description
@ -42,6 +43,22 @@ directions, but it can be limited to the xy-, xz-, yz-plane and the
x-, y-, or z-direction, thus restraining the atoms to a line or a
plane, respectively.
The force constant *k* can be specified as an equal-style or atom-style
:doc:`variable <variable>`. 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 evaluated each time step, and its value(s) will be used as
force constant for the spring force.
Equal-style variables can specify formulas with various mathematical
functions and include :doc:`thermo_style <thermo_style>` command
keywords for the simulation box parameters, time step, and elapsed time.
Thus, it is easy to specify a time-dependent force field.
Atom-style variables can specify the same formulas as equal-style
variables but can also include per-atom values, such as atom
coordinates. Thus, it is easy to specify a spatially-dependent force
field with optional time-dependence as well.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@ -89,7 +106,9 @@ invoked by the :doc:`minimize <minimize>` command.
Restrictions
""""""""""""
none
The KOKKOS version, *fix spring/self/kk* may only be used with a constant
value of K, not a variable.
Related commands
""""""""""""""""

9
doc/src/fix_wall_region.rst
View File

@ -1,8 +1,11 @@
.. index:: fix wall/region
.. index:: fix wall/region/kk
fix wall/region command
=======================
Accelerator Variants: *wall/region/kk*
Syntax
""""""
@ -234,6 +237,12 @@ invoked by the :doc:`minimize <minimize>` command.
minimized), you MUST enable the :doc:`fix_modify <fix_modify>`
*energy* option for this fix.
----------
.. include:: accel_styles.rst
----------
Restrictions
""""""""""""
none

4
doc/src/group.rst
View File

@ -162,7 +162,7 @@ potential energy is above the threshold value :math:`-3.0`.
compute 1 all pe/atom
compute 2 all reduce sum c_1
thermo_style custom step temp pe c_2
run 0
run 0 post no
variable eatom atom "c_1 > -3.0"
group hienergy variable eatom
@ -173,7 +173,7 @@ Note that these lines
compute 2 all reduce sum c_1
thermo_style custom step temp pe c_2
run 0
run 0 post no
are necessary to ensure that the "eatom" variable is current when the
group command invokes it. Because the eatom variable computes the

35
doc/src/pair_bpm_spring.rst
View File

@ -8,7 +8,14 @@ Syntax
.. code-block:: LAMMPS
pair_style bpm/spring
pair_style bpm/spring keyword value ...
* optional keyword = *anharmonic*
.. parsed-literal::
*anharmonic* value = *yes* or *no*
whether forces include the anharmonic term
Examples
""""""""
@ -17,7 +24,8 @@ Examples
pair_style bpm/spring
pair_coeff * * 1.0 1.0 1.0
pair_coeff 1 1 1.0 1.0 1.0
pair_style bpm/spring anharmonic yes
pair_coeff 1 1 1.0 1.0 1.0 50.0
Description
"""""""""""
@ -28,12 +36,16 @@ Style *bpm/spring* computes pairwise forces with the formula
.. math::
F = k (r - r_c)
F = k (r - r_c) + k_a (r - r_c)^3
where :math:`k` is a stiffness and :math:`r_c` is the cutoff length.
An additional damping force is also applied to interacting
particles. The force is proportional to the difference in the
normal velocity of particles
where :math:`k` is a stiffness, :math:`r_c` is the cutoff
length, and :math:`k_a` is an optional anharmonic cubic prefactor
that can be enabled using the *anharmonic* keyword. The anharmonic
term may be useful in scenarios that need to prevent large particle overlap.
An additional damping force is also applied to interacting particles.
The force is proportional to the difference in the normal velocity of
particles
.. math::
@ -73,6 +85,12 @@ commands, or by mixing as described below:
* :math:`r_c` (distance units)
* :math:`\gamma` (force/velocity units)
.. versionadded:: TBD
Additionally, if *anharmonic* is set to *yes*, a fourth coefficient
must be provided:
* :math:`k_a` (force/distance\^3 units)
----------
@ -117,4 +135,5 @@ Related commands
Default
"""""""
none
The option defaults are *anharmonic* = *no*

2
doc/src/pair_coul.rst
View File

@ -180,7 +180,7 @@ coulomb styles in :doc:`hybrid pair styles <pair_hybrid>`.
----------
.. versionadded:: TBD
.. versionadded:: 19Nov2024
Style *coul/ctip* computes the Coulomb interactions as described in
:ref:`Plummer <Plummer1>`. It uses the the damped shifted model as in

32
doc/src/pair_granular.rst
View File

@ -235,12 +235,15 @@ given by:
.. math::
\eta_n = \alpha (m_{eff}k_n)^{1/2}
\eta_n = \alpha (m_{eff}k_{nd})^{1/2}
For normal contact models based on material parameters, :math:`k_n = 4/3Ea`. This
damping model is not compatible with cohesive normal models such as *JKR* or *DMT*.
The parameter :math:`\alpha` is related to the restitution coefficient *e*
according to:
where :math:`k_{nd}` is an effective harmonic stiffness equal to the ratio of
the normal force to the overlap. For example, :math:`k_{nd} = 4/3Ea` for a
Hertz contact model based on material parameters with :math:`a` being
the contact radius of :math:`\sqrt{\delta R}`. For Hooke, :math:`k_{nd}`
is simply the spring constant or :math:`k_{n}`. This damping model is not
compatible with cohesive normal models such as *JKR* or *DMT*. The parameter
:math:`\alpha` is related to the restitution coefficient *e* according to:
.. math::
@ -251,25 +254,26 @@ of the normal contact model parameters should be between 0 and 1, but
no error check is performed on this.
The *coeff_restitution* model is useful when a specific normal coefficient of
restitution :math:`e` is required. In these models, the normal coefficient of
restitution :math:`e` is specified as an input. Following the approach of
:ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke* normal model,
*coeff_restitution* calculates the damping coefficient as:
restitution :math:`e` is required. It operates much like the *Tsuji* model
but, the normal coefficient of restitution :math:`e` is specified as an input
in place of the usual :math:`\eta_{n0}` value in the normal model. Following
the approach of :ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke*
normal model, *coeff_restitution* then calculates the damping coefficient as:
.. math::
\eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
\eta_n = \sqrt{\frac{4m_{eff}k_{nd}}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
where :math:`k_{nd}` is the same stiffness defined in the above *Tsuji* model.
For any other normal model, e.g. the *hertz* and *hertz/material* models, the damping
coefficient is:
.. math::
\eta_n = -2\sqrt{\frac{5}{6}}\frac{\log(e)}{\sqrt{\pi^2+(\log(e))^2}}(R_{eff} \delta_{ij})^{\frac{1}{4}}\sqrt{\frac{3}{2}k_n m_{eff}} ,
\eta_n = -2\sqrt{\frac{5}{6}}\frac{\log(e)}{\sqrt{\pi^2+(\log(e))^2}}\sqrt{\frac{3}{2}k_{nd} m_{eff}} ,
where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since
*coeff_restitution* accounts for the effective mass, effective radius, and
pairwise overlaps (except when used with the *hooke* normal model) when calculating
Since *coeff_restitution* accounts for the effective mass, effective radius,
and pairwise overlaps (except when used with the *hooke* normal model) when calculating
the damping coefficient, it accurately reproduces the specified coefficient of
restitution for both monodisperse and polydisperse particle pairs. This damping
model is not compatible with cohesive normal models such as *JKR* or *DMT*.

26
doc/src/pair_reaxff.rst
View File

@ -20,7 +20,7 @@ Syntax
.. parsed-literal::
keyword = *checkqeq* or *lgvdw* or *safezone* or *mincap* or *minhbonds* or *tabulate* or *list/blocking*
*checkqeq* value = *yes* or *no* = whether or not to require qeq/reaxff or acks2/reaxff fix
*checkqeq* value = *yes* or *no* = whether or not to require one of fix qeq/reaxff, fix acks2/reaxff or fix qtpie/reaxff
*enobonds* value = *yes* or *no* = whether or not to tally energy of atoms with no bonds
*lgvdw* value = *yes* or *no* = whether or not to use a low gradient vdW correction
*safezone* = factor used for array allocation
@ -120,20 +120,22 @@ up that process.
The ReaxFF parameter files provided were created using a charge
equilibration (QEq) model for handling the electrostatic interactions.
Therefore, by default, LAMMPS requires that either the
:doc:`fix qeq/reaxff <fix_qeq_reaxff>` or the
:doc:`fix qeq/shielded <fix_qeq>` or :doc:`fix acks2/reaxff <fix_acks2_reaxff>`
command be used with
*pair_style reaxff* when simulating a ReaxFF model, to equilibrate
the charges each timestep.
Therefore, by default, LAMMPS requires that
:doc:`fix qeq/reaxff <fix_qeq_reaxff>` or :doc:`fix qeq/shielded <fix_qeq>`
or :doc:`fix acks2/reaxff <fix_acks2_reaxff>`
or :doc:`fix qtpie/reaxff <fix_qtpie_reaxff>`
is used with *pair_style reaxff* when simulating a ReaxFF model,
to equilibrate the charges at each timestep.
See the :doc:`fix qeq/reaxff <fix_qeq_reaxff>` or :doc:`fix qeq/shielded <fix_qeq>`
or :doc:`fix acks2/reaxff <fix_acks2_reaxff>`
or :doc:`fix qtpie/reaxff <fix_qtpie_reaxff>`
command documentation for more details.
Using the keyword *checkqeq* with the value *no* turns off the check
for the QEq fixes, allowing a simulation to be run without charge
equilibration. In this case, the static charges you assign to each
atom will be used for computing the electrostatic interactions in
the system. See the :doc:`fix qeq/reaxff <fix_qeq_reaxff>` or
:doc:`fix qeq/shielded <fix_qeq>` or :doc:`fix acks2/reaxff <fix_acks2_reaxff>`
command documentation for more details.
the system.
Using the optional keyword *lgvdw* with the value *yes* turns on the
low-gradient correction of ReaxFF for long-range London Dispersion,
@ -372,8 +374,8 @@ Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`, :doc:`fix qeq/reaxff <fix_qeq_reaxff>`,
:doc:`fix acks2/reaxff <fix_acks2_reaxff>`, :doc:`fix reaxff/bonds <fix_reaxff_bonds>`,
:doc:`fix reaxff/species <fix_reaxff_species>`,
:doc:`fix acks2/reaxff <fix_acks2_reaxff>`, :doc:`fix qtpie/reaxff <fix_qtpie_reaxff>`,
:doc:`fix reaxff/bonds <fix_reaxff_bonds>`, :doc:`fix reaxff/species <fix_reaxff_species>`,
:doc:`compute reaxff/atom <compute_reaxff_atom>`
Default

6
doc/src/pair_sph_heatconduction.rst
View File

@ -30,6 +30,12 @@ The transport model is the diffusion equation for the internal energy.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above.

6
doc/src/pair_sph_idealgas.rst
View File

@ -36,6 +36,12 @@ particles from interpenetrating :ref:`(Monaghan) <ideal-Monoghan>`.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above.

6
doc/src/pair_sph_lj.rst
View File

@ -34,6 +34,12 @@ interpenetrating :ref:`(Monaghan) <Monoghan>`.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above.

6
doc/src/pair_sph_rhosum.rst
View File

@ -29,6 +29,12 @@ SPH particles by kernel function interpolation, every Nstep timesteps.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above.

6
doc/src/pair_sph_taitwater.rst
View File

@ -41,6 +41,12 @@ prevent particles from interpenetrating :ref:`(Monaghan) <Monaghan>`.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above.

6
doc/src/pair_sph_taitwater_morris.rst
View File

@ -37,6 +37,12 @@ This pair style also computes laminar viscosity :ref:`(Morris) <Morris>`.
See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in
LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above.

8
doc/src/region.rst
View File

@ -3,6 +3,8 @@
region command
==============
Accelerator Variants: *block/kk*, *sphere/kk*
Syntax
""""""
@ -74,7 +76,7 @@ Syntax
Rx,Ry,Rz = axis of rotation vector
*open* value = integer from 1-6 corresponding to face index (see below)
* accelerated styles (with same args) = *block/kk*
* accelerated styles (with same args) = *block/kk, sphere/kk*
Examples
""""""""
@ -401,9 +403,9 @@ sub-regions can be defined with the *open* keyword.
.. note::
Currently, only *block* style regions are supported by Kokkos. The
Currently, only *block* and *sphere* style regions are supported by KOKKOS. The
code using the region (such as a fix or compute) must also be
supported by Kokkos or no acceleration will occur.
supported by KOKKOS or no acceleration will occur.
----------

6
doc/src/set.rst
View File

@ -516,6 +516,12 @@ Keywords *sph/e*, *sph/cv*, and *sph/rho* set the energy, heat capacity,
and density of smoothed particle hydrodynamics (SPH) particles. See
`this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in LAMMPS.
.. note::
Please note that the SPH PDF guide file has not been updated for
many years and thus does not reflect the current *syntax* of the
SPH package commands. For that please refer to the LAMMPS manual.
Keyword *smd/mass/density* sets the mass of all selected particles, but
it is only applicable to the Smooth Mach Dynamics package MACHDYN. It
assumes that the particle volume has already been correctly set and

9
doc/utils/check-styles.py
View File

@ -73,8 +73,6 @@ omp = re.compile("(.+)/omp\\s*$")
opt = re.compile("(.+)/opt\\s*$")
removed = re.compile("(.*)Deprecated$")
accel_pattern = re.compile(r"^.. include::\s+accel_styles.rst$")
def require_accel_include(path):
found = False
needs = False
@ -94,6 +92,7 @@ def require_accel_include(path):
if kokkos.match(line): needs = True
if intel.match(line): needs = True
if opt.match(line): needs = True
if path == "src/fix_colvars.rst": needs = False
m = cmd_pattern.match(line)
if m:
if gpu.match(line): needs = True
@ -167,7 +166,9 @@ def check_style(filename, dirname, pattern, styles, name, suffix=False, skip=set
# known undocumented aliases we need to skip
if c in skip: continue
s = c
if suffix: s = add_suffix(styles, c)
if suffix:
s = add_suffix(styles, c)
if s == 'colvars (k)' : continue
if not s in matches:
if not styles[c]['removed']:
print(f"{name} style entry {s} is missing or incomplete in {filename}")
@ -300,7 +301,7 @@ for command_type, entries in index.items():
print("Total number of style index entries:", total_index)
skip_angle = ('sdk')
skip_fix = ('python', 'NEIGH_HISTORY/omp','acks2/reax','qeq/reax','reax/c/bonds','reax/c/species', 'pimd')
skip_fix = ('python', 'NEIGH_HISTORY/omp','acks2/reax','qeq/reax','reax/c/bonds','reax/c/species', 'pimd', 'colvars/kk')
skip_pair = ('meam/c','lj/sf','reax/c','lj/sdk','lj/sdk/coul/long','lj/sdk/coul/msm')
skip_compute = ('pressure/cylinder')

2
doc/utils/requirements.txt
View File

@ -1,4 +1,4 @@
Sphinx >= 5.3.0, <9.0
Sphinx >= 5.3.0, <8.2.0
sphinxcontrib-spelling
sphinxcontrib-jquery
sphinx-design

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

@ -815,6 +815,7 @@ dipoleflag
dir
Direc
dirname
Dirkmann
discoverable
discretization
discretized
@ -976,6 +977,7 @@ elaplong
elastance
Electroneg
electronegative
electronegativities
electronegativity
electroneutral
electroneutrality
@ -1286,12 +1288,14 @@ gcc
gcmc
gdot
GeC
Geesthacht
Geier
gencode
Geocomputing
georg
Georg
Geotechnica
Gergs
germain
Germann
Germano
@ -1304,6 +1308,7 @@ gettimeofday
geturl
gewald
Gezelter
gfile
Gflop
gfortran
ghostneigh
@ -1711,6 +1716,7 @@ Jewett
jgissing
ji
Jiang
Jiahao
Jiao
jik
JIK
@ -2365,6 +2371,7 @@ mui
Mukherjee
Mulders
Müller
Mulliken
mult
multi
multibody
@ -2389,6 +2396,7 @@ Murdick
Murtola
Murty
Muser
Mussenbrock
mutexes
Muto
muVT
@ -2491,6 +2499,7 @@ neel
Neel
Neelov
Negre
Negrut
nelem
Nelement
Nelements
@ -3072,6 +3081,7 @@ qqr
qqrd
Qsb
qtb
qtpie
quadratically
quadrupolar
quadrupole
@ -3107,6 +3117,7 @@ Rafferty
rahman
Rahman
Rajamanickam
Rakhsha
Ralf
Raman
ramped
@ -4219,6 +4230,7 @@ zeeman
Zeeman
Zemer
zenodo
Zentrum
Zepeda
zflag
Zhang

4
examples/bpm/impact/in.bpm.impact.rotational
View File

@ -45,7 +45,7 @@ thermo 100
thermo_modify lost ignore lost/bond ignore
#dump 1 all custom 100 atomDump id radius x y z c_nbond
dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
dump_modify 2 header no
#dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
#dump_modify 2 header no
run 7500

4
examples/bpm/impact/in.bpm.impact.spring
View File

@ -47,7 +47,7 @@ thermo 100
thermo_modify lost ignore lost/bond ignore
#dump 1 all custom 100 atomDump id x y z c_nbond
dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
dump_modify 2 header no
#dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
#dump_modify 2 header no
run 7500

4
examples/bpm/pour/in.bpm.pour
View File

@ -13,9 +13,9 @@ region wall_cyl cylinder z 0.0 0.0 10.0 EDGE EDGE side in
region dropzone cylinder z 0.0 0.0 10.0 40.0 50.0 side in
pair_style gran/hertz/history 1.0 NULL 0.5 NULL 0.1 1
bond_style bpm/rotational
bond_style bpm/rotational break no smooth no
pair_coeff 1 1
bond_coeff 1 1.0 0.2 0.01 0.01 2.0 0.4 0.02 0.02 0.2 0.04 0.002 0.002
bond_coeff 1 1.0 0.2 0.01 0.01 0.0 0.0 0.0 0.0 0.2 0.04 0.002 0.002
compute nbond all nbond/atom
compute tbond all reduce sum c_nbond

1117
examples/bpm/pour/log.29Aug2024.bpm.pour.g++.4 Normal file
View File

File diff suppressed because it is too large Load Diff

1096
examples/bpm/pour/log.4May2022.bpm.pour.g++.4
View File

File diff suppressed because it is too large Load Diff

2
examples/comb/README
View File

@ -20,5 +20,3 @@ Examples:
4. in.comb.Cu2O.elastic: Cu2O crystal, qeq on, minimizes, then calculates
elastic constants
5. in.comb.HfO2: HfO2 polymorphs: Monoclinic HfO2 NVT @ 300K
6. in.comb.CuaS: Metallic Cu and amorphous silica interface, qeq on,
five step NVE run

4147
examples/comb/data.CuaS
View File

File diff suppressed because it is too large Load Diff

5
examples/reaxff/water/gauss_exp.txt Normal file
View File

@ -0,0 +1,5 @@
# Gaussian orbital exponents (required for fix qtpie/reaxff) taken from Table 2.2
# of Chen, J. (2009). Theory and applications of fluctuating-charge models.
# The units of the exponents are 1 / (Bohr radius)^2 .
1 0.2240 # O
2 0.5434 # H

29
examples/reaxff/water/in.water.qtpie Normal file
View File

@ -0,0 +1,29 @@
# QTPIE Water
boundary p p p
units real
atom_style charge
read_data data.water
variable x index 1
variable y index 1
variable z index 1
replicate $x $y $z
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * qeq_ff.water O H
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 4928459 rot yes dist gaussian
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt
fix 2 all nvt temp 300 300 50.0
timestep 0.5
thermo 10
thermo_style custom step temp press density vol
run 20

30
examples/reaxff/water/in.water.qtpie.field Normal file
View File

@ -0,0 +1,30 @@
# QTPIE Water
boundary p p p
units real
atom_style charge
read_data data.water
variable x index 1
variable y index 1
variable z index 1
replicate $x $y $z
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * qeq_ff.water O H
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 4928459 rot yes dist gaussian
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt
fix 2 all nvt temp 300 300 50.0
fix 3 all efield 0.0 0.0 0.05
timestep 0.5
thermo 10
thermo_style custom step temp press density vol
run 20

127
examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 Normal file
View File

@ -0,0 +1,127 @@
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b)
using 1 OpenMP thread(s) per MPI task
# QTPIE Water
boundary p p p
units real
atom_style charge
read_data data.water
Reading data file ...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 1 by 1 MPI processor grid
reading atoms ...
3000 atoms
read_data CPU = 0.056 seconds
variable x index 1
variable y index 1
variable z index 1
replicate $x $y $z
replicate 1 $y $z
replicate 1 1 $z
replicate 1 1 1
Replication is creating a 1x1x1 = 1 times larger system...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 1 by 1 MPI processor grid
3000 atoms
replicate CPU = 0.001 seconds
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * qeq_ff.water O H
WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294)
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 4928459 rot yes dist gaussian
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt
fix 2 all nvt temp 300 300 50.0
fix 3 all efield 0.0 0.0 0.05
timestep 0.5
thermo 10
thermo_style custom step temp press density vol
run 20
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
@Article{Gissinger24,
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
journal = {J. Phys. Chem. B},
year = 2024,
volume = 128,
number = 13,
pages = {3282–-3297}
}
- pair reaxff command: doi:10.1016/j.parco.2011.08.005
@Article{Aktulga12,
author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama},
title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques},
journal = {Parallel Computing},
year = 2012,
volume = 38,
number = {4--5},
pages = {245--259}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.5
ghost atom cutoff = 10.5
binsize = 5.25, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair reaxff, perpetual
attributes: half, newton off, ghost
pair build: half/bin/ghost/newtoff
stencil: full/ghost/bin/3d
bin: standard
(2) fix qtpie/reaxff, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 539.2 | 539.2 | 539.2 Mbytes
Step Temp Press Density Volume
0 300 10137.041 1 29915.273
10 296.09128 3564.7969 1 29915.273
20 293.04308 10299.201 1 29915.273
Loop time of 10.7863 on 1 procs for 20 steps with 3000 atoms
Performance: 0.080 ns/day, 299.620 hours/ns, 1.854 timesteps/s, 5.563 katom-step/s
100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 4.7275 | 4.7275 | 4.7275 | 0.0 | 43.83
Neigh | 0.17533 | 0.17533 | 0.17533 | 0.0 | 1.63
Comm | 0.0017376 | 0.0017376 | 0.0017376 | 0.0 | 0.02
Output | 8.2065e-05 | 8.2065e-05 | 8.2065e-05 | 0.0 | 0.00
Modify | 5.8812 | 5.8812 | 5.8812 | 0.0 | 54.52
Other | | 0.0005226 | | | 0.00
Nlocal: 3000 ave 3000 max 3000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 11077 ave 11077 max 11077 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 971775 ave 971775 max 971775 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 971775
Ave neighs/atom = 323.925
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:12

127
examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 Normal file
View File

@ -0,0 +1,127 @@
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b)
using 1 OpenMP thread(s) per MPI task
# QTPIE Water
boundary p p p
units real
atom_style charge
read_data data.water
Reading data file ...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 2 by 2 MPI processor grid
reading atoms ...
3000 atoms
read_data CPU = 0.053 seconds
variable x index 1
variable y index 1
variable z index 1
replicate $x $y $z
replicate 1 $y $z
replicate 1 1 $z
replicate 1 1 1
Replication is creating a 1x1x1 = 1 times larger system...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 2 by 2 MPI processor grid
3000 atoms
replicate CPU = 0.002 seconds
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * qeq_ff.water O H
WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294)
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 4928459 rot yes dist gaussian
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt
fix 2 all nvt temp 300 300 50.0
fix 3 all efield 0.0 0.0 0.05
timestep 0.5
thermo 10
thermo_style custom step temp press density vol
run 20
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
@Article{Gissinger24,
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
journal = {J. Phys. Chem. B},
year = 2024,
volume = 128,
number = 13,
pages = {3282–-3297}
}
- pair reaxff command: doi:10.1016/j.parco.2011.08.005
@Article{Aktulga12,
author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama},
title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques},
journal = {Parallel Computing},
year = 2012,
volume = 38,
number = {4--5},
pages = {245--259}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.5
ghost atom cutoff = 10.5
binsize = 5.25, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair reaxff, perpetual
attributes: half, newton off, ghost
pair build: half/bin/ghost/newtoff
stencil: full/ghost/bin/3d
bin: standard
(2) fix qtpie/reaxff, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 260.5 | 262.2 | 263.6 Mbytes
Step Temp Press Density Volume
0 300 10137.041 1 29915.273
10 296.09128 3564.7969 1 29915.273
20 293.04308 10299.201 1 29915.273
Loop time of 3.14492 on 4 procs for 20 steps with 3000 atoms
Performance: 0.275 ns/day, 87.359 hours/ns, 6.359 timesteps/s, 19.078 katom-step/s
99.6% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.6557 | 1.6847 | 1.7281 | 2.1 | 53.57
Neigh | 0.086503 | 0.086968 | 0.087627 | 0.2 | 2.77
Comm | 0.003309 | 0.046699 | 0.075729 | 12.4 | 1.48
Output | 5.0156e-05 | 5.483e-05 | 6.8111e-05 | 0.0 | 0.00
Modify | 1.3254 | 1.3261 | 1.3266 | 0.0 | 42.16
Other | | 0.0004552 | | | 0.01
Nlocal: 750 ave 760 max 735 min
Histogram: 1 0 0 0 1 0 0 0 0 2
Nghost: 6230.5 ave 6253 max 6193 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Neighs: 276995 ave 280886 max 271360 min
Histogram: 1 0 0 0 1 0 0 0 1 1
Total # of neighbors = 1107981
Ave neighs/atom = 369.327
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:03

126
examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 Normal file
View File

@ -0,0 +1,126 @@
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b)
using 1 OpenMP thread(s) per MPI task
# QTPIE Water
boundary p p p
units real
atom_style charge
read_data data.water
Reading data file ...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 1 by 1 MPI processor grid
reading atoms ...
3000 atoms
read_data CPU = 0.055 seconds
variable x index 1
variable y index 1
variable z index 1
replicate $x $y $z
replicate 1 $y $z
replicate 1 1 $z
replicate 1 1 1
Replication is creating a 1x1x1 = 1 times larger system...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 1 by 1 MPI processor grid
3000 atoms
replicate CPU = 0.001 seconds
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * qeq_ff.water O H
WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294)
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 4928459 rot yes dist gaussian
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt
fix 2 all nvt temp 300 300 50.0
timestep 0.5
thermo 10
thermo_style custom step temp press density vol
run 20
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
@Article{Gissinger24,
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
journal = {J. Phys. Chem. B},
year = 2024,
volume = 128,
number = 13,
pages = {3282–-3297}
}
- pair reaxff command: doi:10.1016/j.parco.2011.08.005
@Article{Aktulga12,
author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama},
title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques},
journal = {Parallel Computing},
year = 2012,
volume = 38,
number = {4--5},
pages = {245--259}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.5
ghost atom cutoff = 10.5
binsize = 5.25, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair reaxff, perpetual
attributes: half, newton off, ghost
pair build: half/bin/ghost/newtoff
stencil: full/ghost/bin/3d
bin: standard
(2) fix qtpie/reaxff, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 539.2 | 539.2 | 539.2 Mbytes
Step Temp Press Density Volume
0 300 10138.375 1 29915.273
10 295.97879 3575.2769 1 29915.273
20 292.76583 10309.128 1 29915.273
Loop time of 10.8138 on 1 procs for 20 steps with 3000 atoms
Performance: 0.080 ns/day, 300.383 hours/ns, 1.849 timesteps/s, 5.548 katom-step/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 4.7177 | 4.7177 | 4.7177 | 0.0 | 43.63
Neigh | 0.17607 | 0.17607 | 0.17607 | 0.0 | 1.63
Comm | 0.0017295 | 0.0017295 | 0.0017295 | 0.0 | 0.02
Output | 8.5431e-05 | 8.5431e-05 | 8.5431e-05 | 0.0 | 0.00
Modify | 5.9177 | 5.9177 | 5.9177 | 0.0 | 54.72
Other | | 0.0004911 | | | 0.00
Nlocal: 3000 ave 3000 max 3000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 11077 ave 11077 max 11077 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 971830 ave 971830 max 971830 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 971830
Ave neighs/atom = 323.94333
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:12

126
examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 Normal file
View File

@ -0,0 +1,126 @@
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b)
using 1 OpenMP thread(s) per MPI task
# QTPIE Water
boundary p p p
units real
atom_style charge
read_data data.water
Reading data file ...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 2 by 2 MPI processor grid
reading atoms ...
3000 atoms
read_data CPU = 0.053 seconds
variable x index 1
variable y index 1
variable z index 1
replicate $x $y $z
replicate 1 $y $z
replicate 1 1 $z
replicate 1 1 1
Replication is creating a 1x1x1 = 1 times larger system...
orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046)
1 by 2 by 2 MPI processor grid
3000 atoms
replicate CPU = 0.002 seconds
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * qeq_ff.water O H
WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294)
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
velocity all create 300.0 4928459 rot yes dist gaussian
fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt
fix 2 all nvt temp 300 300 50.0
timestep 0.5
thermo 10
thermo_style custom step temp press density vol
run 20
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
@Article{Gissinger24,
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
journal = {J. Phys. Chem. B},
year = 2024,
volume = 128,
number = 13,
pages = {3282–-3297}
}
- pair reaxff command: doi:10.1016/j.parco.2011.08.005
@Article{Aktulga12,
author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama},
title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques},
journal = {Parallel Computing},
year = 2012,
volume = 38,
number = {4--5},
pages = {245--259}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.5
ghost atom cutoff = 10.5
binsize = 5.25, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair reaxff, perpetual
attributes: half, newton off, ghost
pair build: half/bin/ghost/newtoff
stencil: full/ghost/bin/3d
bin: standard
(2) fix qtpie/reaxff, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 260.5 | 262.2 | 263.6 Mbytes
Step Temp Press Density Volume
0 300 10138.375 1 29915.273
10 295.97879 3575.2769 1 29915.273
20 292.76583 10309.128 1 29915.273
Loop time of 3.13598 on 4 procs for 20 steps with 3000 atoms
Performance: 0.276 ns/day, 87.111 hours/ns, 6.378 timesteps/s, 19.133 katom-step/s
99.6% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.6622 | 1.695 | 1.7252 | 2.2 | 54.05
Neigh | 0.086543 | 0.087117 | 0.087848 | 0.2 | 2.78
Comm | 0.0048192 | 0.035002 | 0.067754 | 15.4 | 1.12
Output | 4.8033e-05 | 5.3375e-05 | 6.6893e-05 | 0.0 | 0.00
Modify | 1.3176 | 1.3183 | 1.3189 | 0.0 | 42.04
Other | | 0.0004753 | | | 0.02
Nlocal: 750 ave 760 max 735 min
Histogram: 1 0 0 0 1 0 0 0 0 2
Nghost: 6229.5 ave 6253 max 6191 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Neighs: 277011 ave 280900 max 271380 min
Histogram: 1 0 0 0 1 0 0 0 1 1
Total # of neighbors = 1108044
Ave neighs/atom = 369.348
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:03

102
examples/threebody/log.08Oct24.mos2.sw.mod.g++.4 Normal file
View File

@ -0,0 +1,102 @@
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-512-g13c57ab9b5)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# monolayer MoS2
units metal
boundary p p f
processors * * 1
atom_style atomic
read_data single_layer_MoS2.data
Reading data file ...
triclinic box = (0 0 -100) to (51.15232 44.299209 100) with tilt (25.57616 0 0)
WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (src/domain.cpp:221)
2 by 2 by 1 MPI processor grid
reading atoms ...
768 atoms
read_data CPU = 0.003 seconds
mass * 32.065 # mass of sulphur atom , uint: a.u.=1.66X10^(-27)kg
mass 1 95.94 # mass of molebdenum atom , uint: a.u.=1.66X10^(-27)kg
########################## Define potentials ################################
pair_style sw/mod maxdelcs 0.25 0.35
pair_coeff * * tmd.sw.mod Mo S S
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
#########################################################################
### Simulation settings ####
timestep 0.001
velocity all create 300.0 12345
############################
# Output
thermo 500
thermo_style custom step etotal pe ke temp
thermo_modify lost warn
###### Run molecular dynamics ######
fix thermostat all nve
run 5000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.158796
ghost atom cutoff = 5.158796
binsize = 2.579398, bins = 30 18 78
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair sw/mod, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.165 | 3.165 | 3.165 Mbytes
Step TotEng PotEng KinEng Temp
0 -899.28605 -929.02881 29.742759 300
500 -899.28626 -922.45519 23.168929 233.69313
1000 -899.29247 -925.86547 26.573002 268.02828
1500 -899.27957 -916.95478 17.675214 178.28084
2000 -899.28171 -918.38728 19.105573 192.70814
2500 -899.28732 -922.50423 23.21691 234.17709
3000 -899.28195 -918.74112 19.459174 196.27473
3500 -899.27944 -918.03105 18.751604 189.13784
4000 -899.28397 -920.50737 21.223397 214.06955
4500 -899.28386 -919.79154 20.507685 206.85053
5000 -899.28077 -918.78947 19.508698 196.77425
Loop time of 0.595509 on 4 procs for 5000 steps with 768 atoms
Performance: 725.430 ns/day, 0.033 hours/ns, 8396.182 timesteps/s, 6.448 Matom-step/s
99.9% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.4603 | 0.49732 | 0.54269 | 4.2 | 83.51
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.03293 | 0.078347 | 0.11558 | 10.6 | 13.16
Output | 0.00010079 | 0.00010935 | 0.00012827 | 0.0 | 0.02
Modify | 0.0073413 | 0.0082665 | 0.0091767 | 0.7 | 1.39
Other | | 0.01146 | | | 1.92
Nlocal: 192 ave 194 max 190 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Nghost: 194 ave 196 max 192 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 5120 ave 5170 max 5070 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Total # of neighbors = 20480
Ave neighs/atom = 26.666667
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00

28935
examples/tracker/contact_history.dat
View File

File diff suppressed because it is too large Load Diff

7
examples/tracker/in.track
View File

@ -11,7 +11,7 @@ mass 1 1.0
mass 2 1.0
velocity all create 2.0 87287 loop geom
pair_style hybrid/overlay lj/cut 2.5 tracker
pair_style hybrid/overlay lj/cut 2.5 tracker trackfix 1000 id1 id2 time/created time/broken time/total time/min 0.5
pair_coeff * * lj/cut 1.0 1.0 2.5
pair_coeff * * tracker 2.5
@ -19,10 +19,9 @@ neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
fix 2 all pair/tracker 1000 id1 id2 time/created time/broken time/total time/min 0.5
fix 3 all ave/histo 1000 1 1000 0 6 30 f_2[5] mode vector kind local file lifetime_hist.dat
fix 2 all ave/histo 1000 1 1000 0 6 30 f_trackfix[5] mode vector kind local file lifetime_hist.dat
dump 1 all local 1000 contact_history.dat f_2[1] f_2[2] f_2[3] f_2[4]
dump 1 all local 1000 contact_history.dat f_trackfix[1] f_trackfix[2] f_trackfix[3] f_trackfix[4]
dump_modify 1 header no
# compute 1 all property/local patom1 patom2

54
examples/tracker/lifetime_hist.dat
View File

@ -1,4 +1,4 @@
# Histogrammed data for fix 3
# Histogrammed data for fix 2
# TimeStep Number-of-bins Total-counts Missing-counts Min-value Max-value
# Bin Coord Count Count/Total
0 30 0 0 1e+20 -1e+20
@ -32,34 +32,34 @@
28 5.5 0 0
29 5.7 0 0
30 5.9 0 0
1000 30 8122 0 0.5 5
1000 30 528 20285 1 1000
1 0.1 0 0
2 0.3 0 0
3 0.5 910 0.112041
4 0.7 1253 0.154272
5 0.9 953 0.117336
6 1.1 747 0.0919724
7 1.3 559 0.0688254
8 1.5 501 0.0616843
9 1.7 421 0.0518345
10 1.9 356 0.0438316
11 2.1 300 0.0369367
12 2.3 281 0.0345974
13 2.5 242 0.0297956
14 2.7 226 0.0278257
15 2.9 175 0.0215464
16 3.1 168 0.0206846
17 3.3 162 0.0199458
18 3.5 129 0.0158828
19 3.7 151 0.0185915
20 3.9 137 0.0168678
21 4.1 98 0.012066
22 4.3 104 0.0128047
23 4.5 83 0.0102192
24 4.7 77 0.00948042
25 4.9 86 0.0105885
26 5.1 3 0.000369367
3 0.5 0 0
4 0.7 0 0
5 0.9 0 0
6 1.1 17 0.032197
7 1.3 0 0
8 1.5 0 0
9 1.7 0 0
10 1.9 0 0
11 2.1 51 0.0965909
12 2.3 0 0
13 2.5 0 0
14 2.7 0 0
15 2.9 0 0
16 3.1 80 0.151515
17 3.3 0 0
18 3.5 0 0
19 3.7 0 0
20 3.9 0 0
21 4.1 91 0.172348
22 4.3 0 0
23 4.5 0 0
24 4.7 0 0
25 4.9 0 0
26 5.1 117 0.221591
27 5.3 0 0
28 5.5 0 0
29 5.7 0 0
30 5.9 0 0
30 5.9 172 0.325758

108
examples/tracker/log.29Aug24.tracker.g++.1 Normal file
View File

@ -0,0 +1,108 @@
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-731-g5ce635757f)
# 3d Lennard-Jones melt with tracking
units lj
atom_style atomic
lattice fcc 0.8442
Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962
region box block 0 4 0 4 0 4
create_box 2 box
Created orthogonal box = (0 0 0) to (6.7183848 6.7183848 6.7183848)
1 by 1 by 1 MPI processor grid
create_atoms 2 box
Created 256 atoms
using lattice units in orthogonal box = (0 0 0) to (6.7183848 6.7183848 6.7183848)
create_atoms CPU = 0.000 seconds
mass 1 1.0
mass 2 1.0
velocity all create 2.0 87287 loop geom
pair_style hybrid/overlay lj/cut 2.5 tracker trackfix 1000 id1 id2 time/created time/broken time/total time/min 0.5
pair_coeff * * lj/cut 1.0 1.0 2.5
pair_coeff * * tracker 2.5
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
fix 2 all ave/histo 1000 1 1000 0 6 30 f_trackfix[5] mode vector kind local file lifetime_hist.dat
dump 1 all local 1000 contact_history.dat f_trackfix[1] f_trackfix[2] f_trackfix[3] f_trackfix[4]
dump_modify 1 header no
# compute 1 all property/local patom1 patom2
# dump 2 all local 1 pairs.dat c_1[1] c_1[2]
# dump 3 all atom 1 atoms.dat
thermo 50
run 1000 #0
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 20 steps, delay = 0 steps, check = no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
(2) pair tracker, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 8.679 | 8.679 | 8.679 Mbytes
Step Temp E_pair E_mol TotEng Press
0 2 -6.7733681 0 -3.7850868 -4.5535126
50 1.0540428 -5.3670859 0 -3.7921977 2.3386375
100 1.0402713 -5.3493439 0 -3.7950323 2.4553748
150 1.0570745 -5.3738961 0 -3.7944781 2.3767396
200 1.0431846 -5.3518647 0 -3.7932002 2.5010135
250 1.070121 -5.3902744 0 -3.791363 2.4908658
300 1.0667723 -5.3866302 0 -3.7927224 2.3589344
350 1.000601 -5.2859643 0 -3.7909257 2.9065274
400 0.99256113 -5.2738812 0 -3.7908553 2.8595867
450 1.0482542 -5.357452 0 -3.7912128 2.4707397
500 1.0196176 -5.3123538 0 -3.7889017 2.7230338
550 0.98274535 -5.2586303 0 -3.7902706 2.9156947
600 1.0683914 -5.3863229 0 -3.789996 2.3002719
650 1.0130779 -5.303917 0 -3.7902362 2.8726423
700 1.0583333 -5.3737358 0 -3.792437 2.5770307
750 0.98274506 -5.2612464 0 -3.792887 2.9447027
800 1.0294191 -5.332001 0 -3.7939042 2.5293193
850 0.99240027 -5.2735754 0 -3.7907899 2.7672711
900 1.0293488 -5.3306241 0 -3.7926323 2.6054041
950 0.97137182 -5.2424403 0 -3.7910742 3.129989
1000 1.0009431 -5.2864286 0 -3.7908788 2.7536598
Loop time of 0.286151 on 1 procs for 1000 steps with 256 atoms
Performance: 1509690.552 tau/day, 3494.654 timesteps/s, 894.631 katom-step/s
99.4% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.21578 | 0.21578 | 0.21578 | 0.0 | 75.41
Neigh | 0.020262 | 0.020262 | 0.020262 | 0.0 | 7.08
Comm | 0.0051679 | 0.0051679 | 0.0051679 | 0.0 | 1.81
Output | 0.021204 | 0.021204 | 0.021204 | 0.0 | 7.41
Modify | 0.02307 | 0.02307 | 0.02307 | 0.0 | 8.06
Other | | 0.0006654 | | | 0.23
Nlocal: 256 ave 256 max 256 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1327 ave 1327 max 1327 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 9612 ave 9612 max 9612 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 9612
Ave neighs/atom = 37.546875
Neighbor list builds = 50
Dangerous builds not checked
Total wall time: 0:00:00

107
examples/tracker/log.8Apr21.tracker.g++.1
View File

@ -1,107 +0,0 @@
LAMMPS (8 Apr 2021)
# 3d Lennard-Jones melt with tracking
units lj
atom_style atomic
lattice fcc 0.8442
Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962
region box block 0 4 0 4 0 4
create_box 2 box
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (6.7183848 6.7183848 6.7183848)
1 by 1 by 1 MPI processor grid
create_atoms 2 box
Created 256 atoms
create_atoms CPU = 0.001 seconds
mass 1 1.0
mass 2 1.0
velocity all create 2.0 87287 loop geom
pair_style hybrid/overlay lj/cut 2.5 tracker
pair_coeff * * lj/cut 1.0 1.0 2.5
pair_coeff * * tracker 2.5
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no
fix 1 all nve
fix 2 all pair/tracker 1000 id1 id2 time/created time/broken time/total time/min 0.5
fix 3 all ave/histo 1000 1 1000 0 6 30 f_2[5] mode vector kind local file lifetime_hist.dat
dump 1 all local 1000 contact_history.dat f_2[1] f_2[2] f_2[3] f_2[4]
dump_modify 1 header no
# compute 1 all property/local patom1 patom2
# dump 2 all local 1 pairs.dat c_1[1] c_1[2]
# dump 3 all atom 1 atoms.dat
thermo 50
run 1000 #0
Neighbor list info ...
update every 20 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 2.8
ghost atom cutoff = 2.8
binsize = 1.4, bins = 5 5 5
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
(2) pair tracker, perpetual
attributes: half, newton on, history
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 10.46 | 10.46 | 10.46 Mbytes
Step Temp E_pair E_mol TotEng Press
0 2 -6.7733681 0 -3.7850868 -4.5535126
50 1.0540428 -5.3670859 0 -3.7921977 2.3386375
100 1.0402713 -5.3493439 0 -3.7950323 2.4553748
150 1.0570745 -5.3738961 0 -3.7944781 2.3767396
200 1.0431846 -5.3518647 0 -3.7932002 2.5010135
250 1.070121 -5.3902744 0 -3.791363 2.4908658
300 1.0667723 -5.3866302 0 -3.7927224 2.3589344
350 1.000601 -5.2859643 0 -3.7909257 2.9065274
400 0.99256113 -5.2738812 0 -3.7908553 2.8595867
450 1.0482542 -5.357452 0 -3.7912128 2.4707397
500 1.0196176 -5.3123538 0 -3.7889017 2.7230338
550 0.98274535 -5.2586303 0 -3.7902706 2.9156947
600 1.0683914 -5.3863229 0 -3.789996 2.3002719
650 1.0130779 -5.303917 0 -3.7902362 2.8726423
700 1.0583333 -5.3737358 0 -3.792437 2.5770307
750 0.98274506 -5.2612464 0 -3.792887 2.9447027
800 1.0294191 -5.332001 0 -3.7939042 2.5293193
850 0.99240027 -5.2735754 0 -3.7907899 2.7672711
900 1.0293488 -5.3306241 0 -3.7926323 2.6054041
950 0.97137182 -5.2424403 0 -3.7910742 3.129989
1000 1.0009431 -5.2864286 0 -3.7908788 2.7536598
Loop time of 0.310363 on 1 procs for 1000 steps with 256 atoms
Performance: 1391918.252 tau/day, 3222.033 timesteps/s
100.7% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.22132 | 0.22132 | 0.22132 | 0.0 | 71.31
Neigh | 0.03458 | 0.03458 | 0.03458 | 0.0 | 11.14
Comm | 0.0087938 | 0.0087938 | 0.0087938 | 0.0 | 2.83
Output | 0.014075 | 0.014075 | 0.014075 | 0.0 | 4.54
Modify | 0.02818 | 0.02818 | 0.02818 | 0.0 | 9.08
Other | | 0.003411 | | | 1.10
Nlocal: 256.000 ave 256 max 256 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1327.00 ave 1327 max 1327 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 9612.00 ave 9612 max 9612 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 9612
Ave neighs/atom = 37.546875
Neighbor list builds = 50
Dangerous builds not checked
Total wall time: 0:00:00

40
examples/wall/in.wall.block Normal file
View File

@ -0,0 +1,40 @@
units real
molecule water tip3p.mol
atom_style full
variable radius equal 100.0
region box block $(-v_radius) $(v_radius) $(-v_radius) $(v_radius) $(-v_radius) $(v_radius)
create_box 2 box bond/types 1 angle/types 1 &
extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
mass 1 15.9994
mass 2 1.008
bond_style zero
bond_coeff 1 0.9574
angle_style zero
angle_coeff 1 104.52
region block1 block -80 80 -80 80 -80 80
region block2 block -70 70 -70 70 -70 70
create_atoms 0 random 5000 12345 block2 mol water 12345 overlap 2
thermo 1
thermo_style custom step time spcpu temp press etotal pe
fix wall all wall/region block1 harmonic 1000.0 0.0 2.5
fix_modify wall energy yes
pair_style lj/cut/coul/cut 8.0
pair_coeff 1 1 0.1521 3.1507
pair_coeff 2 2 0.0 1.0
velocity all create 300.0 12345
fix 1 all nvt temp 300 300 100.0
fix 2 all shake 0.001 10 10000 b 1 a 1
dump 2 all movie 10 wall.block.mpg type type size 1500 1500 fsaa yes
dump_modify 2 pad 4 acolor * white/red/green/blue/aqua/magenta
run 100

40
examples/wall/in.wall.sphere Normal file
View File

@ -0,0 +1,40 @@
units real
molecule water tip3p.mol
atom_style full
variable radius equal 100.0
region box block $(-v_radius) $(v_radius) $(-v_radius) $(v_radius) $(-v_radius) $(v_radius)
create_box 2 box bond/types 1 angle/types 1 &
extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
mass 1 15.9994
mass 2 1.008
bond_style zero
bond_coeff 1 0.9574
angle_style zero
angle_coeff 1 104.52
region sphere1 sphere 0 0 0 $(v_radius-10) side in
region sphere2 sphere 0 0 0 $(v_radius-20) side in
create_atoms 0 random 5000 12345 sphere2 mol water 12345 overlap 2
thermo 1
thermo_style custom step time spcpu temp press etotal pe
fix wall all wall/region sphere1 harmonic 1000.0 0.0 2.5
fix_modify wall energy yes
pair_style lj/cut/coul/cut 8.0
pair_coeff 1 1 0.1521 3.1507
pair_coeff 2 2 0.0 1.0
velocity all create 300.0 12345
fix 1 all nvt temp 300 300 100.0
fix 2 all shake 0.001 10 10000 b 1 a 1
dump 2 all movie 10 wall.sphere.mpg type type size 1500 1500 fsaa yes
dump_modify 2 pad 4 acolor * white/red/green/blue/aqua/magenta
run 100

62
examples/wall/tip3p.mol Normal file
View File

@ -0,0 +1,62 @@
# Water molecule. TIP3P geometry
3 atoms
2 bonds
1 angles
Coords
1 0.00000 -0.06556 0.00000
2 0.75695 0.52032 0.00000
3 -0.75695 0.52032 0.00000
Types
1 1
2 2
3 2
Charges
1 -0.834
2 0.417
3 0.417
Bonds
1 1 1 2
2 1 1 3
Angles
1 1 2 1 3
Shake Flags
1 1
2 1
3 1
Shake Atoms
1 1 2 3
2 1 2 3
3 1 2 3
Shake Bond Types
1 1 1 1
2 1 1 1
3 1 1 1
Special Bond Counts
1 2 0 0
2 1 1 0
3 1 1 0
Special Bonds
1 2 3
2 1 3
3 1 2

4
lib/gpu/lal_neighbor.cpp
View File

@ -365,7 +365,9 @@ void Neighbor::get_host(const int inum, int *ilist, int *numj,
int i=ilist[ii];
three_ilist[i] = ii;
}
three_ilist.update_device(inum,true);
// needs to transfer _max_atoms because three_ilist indexes all the atoms (local and ghost)
// not just inum (number of neighbor list items)
three_ilist.update_device(_max_atoms,true);
}
time_nbor.stop();

282
lib/linalg/dbdsdc.cpp Normal file
View File

@ -0,0 +1,282 @@
#ifdef __cplusplus
extern "C" {
#endif
#include "lmp_f2c.h"
static integer c__9 = 9;
static integer c__0 = 0;
static doublereal c_b15 = 1.;
static integer c__1 = 1;
static doublereal c_b29 = 0.;
int dbdsdc_(char *uplo, char *compq, integer *n, doublereal *d__, doublereal *e, doublereal *u,
integer *ldu, doublereal *vt, integer *ldvt, doublereal *q, integer *iq,
doublereal *work, integer *iwork, integer *info, ftnlen uplo_len, ftnlen compq_len)
{
integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
doublereal d__1;
double d_lmp_sign(doublereal *, doublereal *), log(doublereal);
integer i__, j, k;
doublereal p, r__;
integer z__, ic, ii, kk;
doublereal cs;
integer is, iu;
doublereal sn;
integer nm1;
doublereal eps;
integer ivt, difl, difr, ierr, perm, mlvl, sqre;
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern int dlasr_(char *, char *, char *, integer *, integer *, doublereal *, doublereal *,
doublereal *, integer *, ftnlen, ftnlen, ftnlen),
dcopy_(integer *, doublereal *, integer *, doublereal *, integer *),
dswap_(integer *, doublereal *, integer *, doublereal *, integer *);
integer poles, iuplo, nsize, start;
extern int dlasd0_(integer *, integer *, doublereal *, doublereal *, doublereal *, integer *,
doublereal *, integer *, integer *, integer *, doublereal *, integer *);
extern doublereal dlamch_(char *, ftnlen);
extern int dlasda_(integer *, integer *, integer *, integer *, doublereal *, doublereal *,
doublereal *, integer *, doublereal *, integer *, doublereal *, doublereal *,
doublereal *, doublereal *, integer *, integer *, integer *, integer *,
doublereal *, doublereal *, doublereal *, doublereal *, integer *,
integer *),
dlascl_(char *, integer *, integer *, doublereal *, doublereal *, integer *, integer *,
doublereal *, integer *, integer *, ftnlen),
dlasdq_(char *, integer *, integer *, integer *, integer *, integer *, doublereal *,
doublereal *, doublereal *, integer *, doublereal *, integer *, doublereal *,
integer *, doublereal *, integer *, ftnlen),
dlaset_(char *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *,
ftnlen),
dlartg_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *,
ftnlen, ftnlen);
extern int xerbla_(char *, integer *, ftnlen);
integer givcol;
extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *, ftnlen);
integer icompq;
doublereal orgnrm;
integer givnum, givptr, qstart, smlsiz, wstart, smlszp;
--d__;
--e;
u_dim1 = *ldu;
u_offset = 1 + u_dim1;
u -= u_offset;
vt_dim1 = *ldvt;
vt_offset = 1 + vt_dim1;
vt -= vt_offset;
--q;
--iq;
--work;
--iwork;
*info = 0;
iuplo = 0;
if (lsame_(uplo, (char *)"U", (ftnlen)1, (ftnlen)1)) {
iuplo = 1;
}
if (lsame_(uplo, (char *)"L", (ftnlen)1, (ftnlen)1)) {
iuplo = 2;
}
if (lsame_(compq, (char *)"N", (ftnlen)1, (ftnlen)1)) {
icompq = 0;
} else if (lsame_(compq, (char *)"P", (ftnlen)1, (ftnlen)1)) {
icompq = 1;
} else if (lsame_(compq, (char *)"I", (ftnlen)1, (ftnlen)1)) {
icompq = 2;
} else {
icompq = -1;
}
if (iuplo == 0) {
*info = -1;
} else if (icompq < 0) {
*info = -2;
} else if (*n < 0) {
*info = -3;
} else if (*ldu < 1 || icompq == 2 && *ldu < *n) {
*info = -7;
} else if (*ldvt < 1 || icompq == 2 && *ldvt < *n) {
*info = -9;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_((char *)"DBDSDC", &i__1, (ftnlen)6);
return 0;
}
if (*n == 0) {
return 0;
}
smlsiz = ilaenv_(&c__9, (char *)"DBDSDC", (char *)" ", &c__0, &c__0, &c__0, &c__0, (ftnlen)6, (ftnlen)1);
if (*n == 1) {
if (icompq == 1) {
q[1] = d_lmp_sign(&c_b15, &d__[1]);
q[smlsiz * *n + 1] = 1.;
} else if (icompq == 2) {
u[u_dim1 + 1] = d_lmp_sign(&c_b15, &d__[1]);
vt[vt_dim1 + 1] = 1.;
}
d__[1] = abs(d__[1]);
return 0;
}
nm1 = *n - 1;
wstart = 1;
qstart = 3;
if (icompq == 1) {
dcopy_(n, &d__[1], &c__1, &q[1], &c__1);
i__1 = *n - 1;
dcopy_(&i__1, &e[1], &c__1, &q[*n + 1], &c__1);
}
if (iuplo == 2) {
qstart = 5;
if (icompq == 2) {
wstart = (*n << 1) - 1;
}
i__1 = *n - 1;
for (i__ = 1; i__ <= i__1; ++i__) {
dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
d__[i__] = r__;
e[i__] = sn * d__[i__ + 1];
d__[i__ + 1] = cs * d__[i__ + 1];
if (icompq == 1) {
q[i__ + (*n << 1)] = cs;
q[i__ + *n * 3] = sn;
} else if (icompq == 2) {
work[i__] = cs;
work[nm1 + i__] = -sn;
}
}
}
if (icompq == 0) {
dlasdq_((char *)"U", &c__0, n, &c__0, &c__0, &c__0, &d__[1], &e[1], &vt[vt_offset], ldvt,
&u[u_offset], ldu, &u[u_offset], ldu, &work[1], info, (ftnlen)1);
goto L40;
}
if (*n <= smlsiz) {
if (icompq == 2) {
dlaset_((char *)"A", n, n, &c_b29, &c_b15, &u[u_offset], ldu, (ftnlen)1);
dlaset_((char *)"A", n, n, &c_b29, &c_b15, &vt[vt_offset], ldvt, (ftnlen)1);
dlasdq_((char *)"U", &c__0, n, n, n, &c__0, &d__[1], &e[1], &vt[vt_offset], ldvt, &u[u_offset],
ldu, &u[u_offset], ldu, &work[wstart], info, (ftnlen)1);
} else if (icompq == 1) {
iu = 1;
ivt = iu + *n;
dlaset_((char *)"A", n, n, &c_b29, &c_b15, &q[iu + (qstart - 1) * *n], n, (ftnlen)1);
dlaset_((char *)"A", n, n, &c_b29, &c_b15, &q[ivt + (qstart - 1) * *n], n, (ftnlen)1);
dlasdq_((char *)"U", &c__0, n, n, n, &c__0, &d__[1], &e[1], &q[ivt + (qstart - 1) * *n], n,
&q[iu + (qstart - 1) * *n], n, &q[iu + (qstart - 1) * *n], n, &work[wstart],
info, (ftnlen)1);
}
goto L40;
}
if (icompq == 2) {
dlaset_((char *)"A", n, n, &c_b29, &c_b15, &u[u_offset], ldu, (ftnlen)1);
dlaset_((char *)"A", n, n, &c_b29, &c_b15, &vt[vt_offset], ldvt, (ftnlen)1);
}
orgnrm = dlanst_((char *)"M", n, &d__[1], &e[1], (ftnlen)1);
if (orgnrm == 0.) {
return 0;
}
dlascl_((char *)"G", &c__0, &c__0, &orgnrm, &c_b15, n, &c__1, &d__[1], n, &ierr, (ftnlen)1);
dlascl_((char *)"G", &c__0, &c__0, &orgnrm, &c_b15, &nm1, &c__1, &e[1], &nm1, &ierr, (ftnlen)1);
eps = dlamch_((char *)"Epsilon", (ftnlen)7) * .9;
mlvl = (integer)(log((doublereal)(*n) / (doublereal)(smlsiz + 1)) / log(2.)) + 1;
smlszp = smlsiz + 1;
if (icompq == 1) {
iu = 1;
ivt = smlsiz + 1;
difl = ivt + smlszp;
difr = difl + mlvl;
z__ = difr + (mlvl << 1);
ic = z__ + mlvl;
is = ic + 1;
poles = is + 1;
givnum = poles + (mlvl << 1);
k = 1;
givptr = 2;
perm = 3;
givcol = perm + mlvl;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
if ((d__1 = d__[i__], abs(d__1)) < eps) {
d__[i__] = d_lmp_sign(&eps, &d__[i__]);
}
}
start = 1;
sqre = 0;
i__1 = nm1;
for (i__ = 1; i__ <= i__1; ++i__) {
if ((d__1 = e[i__], abs(d__1)) < eps || i__ == nm1) {
if (i__ < nm1) {
nsize = i__ - start + 1;
} else if ((d__1 = e[i__], abs(d__1)) >= eps) {
nsize = *n - start + 1;
} else {
nsize = i__ - start + 1;
if (icompq == 2) {
u[*n + *n * u_dim1] = d_lmp_sign(&c_b15, &d__[*n]);
vt[*n + *n * vt_dim1] = 1.;
} else if (icompq == 1) {
q[*n + (qstart - 1) * *n] = d_lmp_sign(&c_b15, &d__[*n]);
q[*n + (smlsiz + qstart - 1) * *n] = 1.;
}
d__[*n] = (d__1 = d__[*n], abs(d__1));
}
if (icompq == 2) {
dlasd0_(&nsize, &sqre, &d__[start], &e[start], &u[start + start * u_dim1], ldu,
&vt[start + start * vt_dim1], ldvt, &smlsiz, &iwork[1], &work[wstart],
info);
} else {
dlasda_(&icompq, &smlsiz, &nsize, &sqre, &d__[start], &e[start],
&q[start + (iu + qstart - 2) * *n], n, &q[start + (ivt + qstart - 2) * *n],
&iq[start + k * *n], &q[start + (difl + qstart - 2) * *n],
&q[start + (difr + qstart - 2) * *n], &q[start + (z__ + qstart - 2) * *n],
&q[start + (poles + qstart - 2) * *n], &iq[start + givptr * *n],
&iq[start + givcol * *n], n, &iq[start + perm * *n],
&q[start + (givnum + qstart - 2) * *n], &q[start + (ic + qstart - 2) * *n],
&q[start + (is + qstart - 2) * *n], &work[wstart], &iwork[1], info);
}
if (*info != 0) {
return 0;
}
start = i__ + 1;
}
}
dlascl_((char *)"G", &c__0, &c__0, &c_b15, &orgnrm, n, &c__1, &d__[1], n, &ierr, (ftnlen)1);
L40:
i__1 = *n;
for (ii = 2; ii <= i__1; ++ii) {
i__ = ii - 1;
kk = i__;
p = d__[i__];
i__2 = *n;
for (j = ii; j <= i__2; ++j) {
if (d__[j] > p) {
kk = j;
p = d__[j];
}
}
if (kk != i__) {
d__[kk] = d__[i__];
d__[i__] = p;
if (icompq == 1) {
iq[i__] = kk;
} else if (icompq == 2) {
dswap_(n, &u[i__ * u_dim1 + 1], &c__1, &u[kk * u_dim1 + 1], &c__1);
dswap_(n, &vt[i__ + vt_dim1], ldvt, &vt[kk + vt_dim1], ldvt);
}
} else if (icompq == 1) {
iq[i__] = i__;
}
}
if (icompq == 1) {
if (iuplo == 1) {
iq[*n] = 1;
} else {
iq[*n] = 0;
}
}
if (iuplo == 2 && icompq == 2) {
dlasr_((char *)"L", (char *)"V", (char *)"B", n, n, &work[1], &work[*n], &u[u_offset], ldu, (ftnlen)1, (ftnlen)1,
(ftnlen)1);
}
return 0;
}
#ifdef __cplusplus
}
#endif

26
lib/linalg/dcombssq.cpp Normal file
View File

@ -0,0 +1,26 @@
#ifdef __cplusplus
extern "C" {
#endif
#include "lmp_f2c.h"
int dcombssq_(doublereal *v1, doublereal *v2)
{
doublereal d__1;
--v2;
--v1;
if (v1[1] >= v2[1]) {
if (v1[1] != 0.) {
d__1 = v2[1] / v1[1];
v1[2] += d__1 * d__1 * v2[2];
} else {
v1[2] += v2[2];
}
} else {
d__1 = v1[1] / v2[1];
v1[2] = v2[2] + d__1 * d__1 * v1[2];
v1[1] = v2[1];
}
return 0;
}
#ifdef __cplusplus
}
#endif

117
lib/linalg/dgebak.cpp Normal file
View File

@ -0,0 +1,117 @@
#ifdef __cplusplus
extern "C" {
#endif
#include "lmp_f2c.h"
int dgebak_(char *job, char *side, integer *n, integer *ilo, integer *ihi, doublereal *scale,
integer *m, doublereal *v, integer *ldv, integer *info, ftnlen job_len, ftnlen side_len)
{
integer v_dim1, v_offset, i__1;
integer i__, k;
doublereal s;
integer ii;
extern int dscal_(integer *, doublereal *, doublereal *, integer *);
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern int dswap_(integer *, doublereal *, integer *, doublereal *, integer *);
logical leftv;
extern int xerbla_(char *, integer *, ftnlen);
logical rightv;
--scale;
v_dim1 = *ldv;
v_offset = 1 + v_dim1;
v -= v_offset;
rightv = lsame_(side, (char *)"R", (ftnlen)1, (ftnlen)1);
leftv = lsame_(side, (char *)"L", (ftnlen)1, (ftnlen)1);
*info = 0;
if (!lsame_(job, (char *)"N", (ftnlen)1, (ftnlen)1) && !lsame_(job, (char *)"P", (ftnlen)1, (ftnlen)1) &&
!lsame_(job, (char *)"S", (ftnlen)1, (ftnlen)1) && !lsame_(job, (char *)"B", (ftnlen)1, (ftnlen)1)) {
*info = -1;
} else if (!rightv && !leftv) {
*info = -2;
} else if (*n < 0) {
*info = -3;
} else if (*ilo < 1 || *ilo > max(1, *n)) {
*info = -4;
} else if (*ihi < min(*ilo, *n) || *ihi > *n) {
*info = -5;
} else if (*m < 0) {
*info = -7;
} else if (*ldv < max(1, *n)) {
*info = -9;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_((char *)"DGEBAK", &i__1, (ftnlen)6);
return 0;
}
if (*n == 0) {
return 0;
}
if (*m == 0) {
return 0;
}
if (lsame_(job, (char *)"N", (ftnlen)1, (ftnlen)1)) {
return 0;
}
if (*ilo == *ihi) {
goto L30;
}
if (lsame_(job, (char *)"S", (ftnlen)1, (ftnlen)1) || lsame_(job, (char *)"B", (ftnlen)1, (ftnlen)1)) {
if (rightv) {
i__1 = *ihi;
for (i__ = *ilo; i__ <= i__1; ++i__) {
s = scale[i__];
dscal_(m, &s, &v[i__ + v_dim1], ldv);
}
}
if (leftv) {
i__1 = *ihi;
for (i__ = *ilo; i__ <= i__1; ++i__) {
s = 1. / scale[i__];
dscal_(m, &s, &v[i__ + v_dim1], ldv);
}
}
}
L30:
if (lsame_(job, (char *)"P", (ftnlen)1, (ftnlen)1) || lsame_(job, (char *)"B", (ftnlen)1, (ftnlen)1)) {
if (rightv) {
i__1 = *n;
for (ii = 1; ii <= i__1; ++ii) {
i__ = ii;
if (i__ >= *ilo && i__ <= *ihi) {
goto L40;
}
if (i__ < *ilo) {
i__ = *ilo - ii;
}
k = (integer)scale[i__];
if (k == i__) {
goto L40;
}
dswap_(m, &v[i__ + v_dim1], ldv, &v[k + v_dim1], ldv);
L40:;
}
}
if (leftv) {
i__1 = *n;
for (ii = 1; ii <= i__1; ++ii) {
i__ = ii;
if (i__ >= *ilo && i__ <= *ihi) {
goto L50;
}
if (i__ < *ilo) {
i__ = *ilo - ii;
}
k = (integer)scale[i__];
if (k == i__) {
goto L50;
}
dswap_(m, &v[i__ + v_dim1], ldv, &v[k + v_dim1], ldv);
L50:;
}
}
}
return 0;
}
#ifdef __cplusplus
}
#endif

513
lib/linalg/dgebal.cpp Normal file
View File

@ -0,0 +1,513 @@
#ifdef __cplusplus
extern "C" {
#endif
#include "lmp_f2c.h"
static integer c__1 = 1;
static integer c__0 = 0;
static integer c_n1 = -1;
int dgebal_(char *job, integer *n, doublereal *a, integer *lda, integer *ilo, integer *ihi,
doublereal *scale, integer *info, ftnlen job_len)
{
integer a_dim1, a_offset, i__1, i__2;
doublereal d__1, d__2;
doublereal c__, f, g;
integer i__, j, k, l, m;
doublereal r__, s, ca, ra;
integer ica, ira, iexc;
extern doublereal dnrm2_(integer *, doublereal *, integer *);
extern int dscal_(integer *, doublereal *, doublereal *, integer *);
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern int dswap_(integer *, doublereal *, integer *, doublereal *, integer *);
doublereal sfmin1, sfmin2, sfmax1, sfmax2;
extern doublereal dlamch_(char *, ftnlen);
extern integer idamax_(integer *, doublereal *, integer *);
extern logical disnan_(doublereal *);
extern int xerbla_(char *, integer *, ftnlen);
logical noconv;
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--scale;
*info = 0;
if (!lsame_(job, (char *)"N", (ftnlen)1, (ftnlen)1) && !lsame_(job, (char *)"P", (ftnlen)1, (ftnlen)1) &&
!lsame_(job, (char *)"S", (ftnlen)1, (ftnlen)1) && !lsame_(job, (char *)"B", (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 *)"DGEBAL", &i__1, (ftnlen)6);
return 0;
}
k = 1;
l = *n;
if (*n == 0) {
goto L210;
}
if (lsame_(job, (char *)"N", (ftnlen)1, (ftnlen)1)) {
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
scale[i__] = 1.;
}
goto L210;
}
if (lsame_(job, (char *)"S", (ftnlen)1, (ftnlen)1)) {
goto L120;
}
goto L50;
L20:
scale[m] = (doublereal)j;
if (j == m) {
goto L30;
}
dswap_(&l, &a[j * a_dim1 + 1], &c__1, &a[m * a_dim1 + 1], &c__1);
i__1 = *n - k + 1;
dswap_(&i__1, &a[j + k * a_dim1], lda, &a[m + k * a_dim1], lda);
L30:
switch (iexc) {
case 1:
goto L40;
case 2:
goto L80;
}
L40:
if (l == 1) {
goto L210;
}
--l;
L50:
for (j = l; j >= 1; --j) {
i__1 = l;
for (i__ = 1; i__ <= i__1; ++i__) {
if (i__ == j) {
goto L60;
}
if (a[j + i__ * a_dim1] != 0.) {
goto L70;
}
L60:;
}
m = l;
iexc = 1;
goto L20;
L70:;
}
goto L90;
L80:
++k;
L90:
i__1 = l;
for (j = k; j <= i__1; ++j) {
i__2 = l;
for (i__ = k; i__ <= i__2; ++i__) {
if (i__ == j) {
goto L100;
}
if (a[i__ + j * a_dim1] != 0.) {
goto L110;
}
L100:;
}
m = k;
iexc = 2;
goto L20;
L110:;
}
L120:
i__1 = l;
for (i__ = k; i__ <= i__1; ++i__) {
scale[i__] = 1.;
}
if (lsame_(job, (char *)"P", (ftnlen)1, (ftnlen)1)) {
goto L210;
}
sfmin1 = dlamch_((char *)"S", (ftnlen)1) / dlamch_((char *)"P", (ftnlen)1);
sfmax1 = 1. / sfmin1;
sfmin2 = sfmin1 * 2.;
sfmax2 = 1. / sfmin2;
L140:
noconv = FALSE_;
i__1 = l;
for (i__ = k; i__ <= i__1; ++i__) {
i__2 = l - k + 1;
c__ = dnrm2_(&i__2, &a[k + i__ * a_dim1], &c__1);
i__2 = l - k + 1;
r__ = dnrm2_(&i__2, &a[i__ + k * a_dim1], lda);
ica = idamax_(&l, &a[i__ * a_dim1 + 1], &c__1);
ca = (d__1 = a[ica + i__ * a_dim1], abs(d__1));
i__2 = *n - k + 1;
ira = idamax_(&i__2, &a[i__ + k * a_dim1], lda);
ra = (d__1 = a[i__ + (ira + k - 1) * a_dim1], abs(d__1));
if (c__ == 0. || r__ == 0.) {
goto L200;
}
g = r__ / 2.;
f = 1.;
s = c__ + r__;
L160:
d__1 = max(f, c__);
d__2 = min(r__, g);
if (c__ >= g || max(d__1, ca) >= sfmax2 || min(d__2, ra) <= sfmin2) {
goto L170;
}
d__1 = c__ + f + ca + r__ + g + ra;
if (disnan_(&d__1)) {
*info = -3;
i__2 = -(*info);
xerbla_((char *)"DGEBAL", &i__2, (ftnlen)6);
return 0;
}
f *= 2.;
c__ *= 2.;
ca *= 2.;
r__ /= 2.;
g /= 2.;
ra /= 2.;
goto L160;
L170:
g = c__ / 2.;
L180:
d__1 = min(f, c__), d__1 = min(d__1, g);
if (g < r__ || max(r__, ra) >= sfmax2 || min(d__1, ca) <= sfmin2) {
goto L190;
}
f /= 2.;
c__ /= 2.;
g /= 2.;
ca /= 2.;
r__ *= 2.;
ra *= 2.;
goto L180;
L190:
if (c__ + r__ >= s * .95) {
goto L200;
}
if (f < 1. && scale[i__] < 1.) {
if (f * scale[i__] <= sfmin1) {
goto L200;
}
}
if (f > 1. && scale[i__] > 1.) {
if (scale[i__] >= sfmax1 / f) {
goto L200;
}
}
g = 1. / f;
scale[i__] *= f;
noconv = TRUE_;
i__2 = *n - k + 1;
dscal_(&i__2, &g, &a[i__ + k * a_dim1], lda);
dscal_(&l, &f, &a[i__ * a_dim1 + 1], &c__1);
L200:;
}
if (noconv) {
goto L140;
}
L210:
*ilo = k;
*ihi = l;
return 0;
}
int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *a, integer *lda, doublereal *wr,
doublereal *wi, doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr,
doublereal *work, integer *lwork, integer *info, ftnlen jobvl_len, ftnlen jobvr_len)
{
integer a_dim1, a_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, i__2, i__3;
doublereal d__1, d__2;
double sqrt(doublereal);
integer i__, k;
doublereal r__, cs, sn;
integer ihi;
doublereal scl;
integer ilo;
doublereal dum[1], eps;
integer lwork_trevc__, ibal;
char side[1];
doublereal anrm;
integer ierr, itau;
extern int drot_(integer *, doublereal *, integer *, doublereal *, integer *, doublereal *,
doublereal *);
integer iwrk, nout;
extern doublereal dnrm2_(integer *, doublereal *, integer *);
extern int dscal_(integer *, doublereal *, doublereal *, integer *);
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern doublereal dlapy2_(doublereal *, doublereal *);
extern int dlabad_(doublereal *, doublereal *),
dgebak_(char *, char *, integer *, integer *, integer *, doublereal *, integer *,
doublereal *, integer *, integer *, ftnlen, ftnlen),
dgebal_(char *, integer *, doublereal *, integer *, integer *, integer *, doublereal *,
integer *, ftnlen);
logical scalea;
extern doublereal dlamch_(char *, ftnlen);
doublereal cscale;
extern doublereal dlange_(char *, integer *, integer *, doublereal *, integer *, doublereal *,
ftnlen);
extern int dgehrd_(integer *, integer *, integer *, doublereal *, integer *, doublereal *,
doublereal *, integer *, integer *),
dlascl_(char *, integer *, integer *, doublereal *, doublereal *, integer *, integer *,
doublereal *, integer *, integer *, ftnlen);
extern integer idamax_(integer *, doublereal *, integer *);
extern int dlacpy_(char *, integer *, integer *, doublereal *, integer *, doublereal *,
integer *, ftnlen),
dlartg_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *),
xerbla_(char *, integer *, ftnlen);
logical select[1];
extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *,
ftnlen, ftnlen);
doublereal bignum;
extern int dorghr_(integer *, integer *, integer *, doublereal *, integer *, doublereal *,
doublereal *, integer *, integer *),
dhseqr_(char *, char *, integer *, integer *, integer *, doublereal *, integer *,
doublereal *, doublereal *, doublereal *, integer *, doublereal *, integer *,
integer *, ftnlen, ftnlen);
integer minwrk, maxwrk;
logical wantvl;
doublereal smlnum;
integer hswork;
logical lquery, wantvr;
extern int dtrevc3_(char *, char *, logical *, integer *, doublereal *, integer *, doublereal *,
integer *, doublereal *, integer *, integer *, integer *, doublereal *,
integer *, integer *, ftnlen, ftnlen);
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--wr;
--wi;
vl_dim1 = *ldvl;
vl_offset = 1 + vl_dim1;
vl -= vl_offset;
vr_dim1 = *ldvr;
vr_offset = 1 + vr_dim1;
vr -= vr_offset;
--work;
*info = 0;
lquery = *lwork == -1;
wantvl = lsame_(jobvl, (char *)"V", (ftnlen)1, (ftnlen)1);
wantvr = lsame_(jobvr, (char *)"V", (ftnlen)1, (ftnlen)1);
if (!wantvl && !lsame_(jobvl, (char *)"N", (ftnlen)1, (ftnlen)1)) {
*info = -1;
} else if (!wantvr && !lsame_(jobvr, (char *)"N", (ftnlen)1, (ftnlen)1)) {
*info = -2;
} else if (*n < 0) {
*info = -3;
} else if (*lda < max(1, *n)) {
*info = -5;
} else if (*ldvl < 1 || wantvl && *ldvl < *n) {
*info = -9;
} else if (*ldvr < 1 || wantvr && *ldvr < *n) {
*info = -11;
}
if (*info == 0) {
if (*n == 0) {
minwrk = 1;
maxwrk = 1;
} else {
maxwrk = (*n << 1) +
*n * ilaenv_(&c__1, (char *)"DGEHRD", (char *)" ", n, &c__1, n, &c__0, (ftnlen)6, (ftnlen)1);
if (wantvl) {
minwrk = *n << 2;
i__1 = maxwrk,
i__2 = (*n << 1) + (*n - 1) * ilaenv_(&c__1, (char *)"DORGHR", (char *)" ", n, &c__1, n, &c_n1,
(ftnlen)6, (ftnlen)1);
maxwrk = max(i__1, i__2);
dhseqr_((char *)"S", (char *)"V", n, &c__1, n, &a[a_offset], lda, &wr[1], &wi[1], &vl[vl_offset],
ldvl, &work[1], &c_n1, info, (ftnlen)1, (ftnlen)1);
hswork = (integer)work[1];
i__1 = maxwrk, i__2 = *n + 1, i__1 = max(i__1, i__2), i__2 = *n + hswork;
maxwrk = max(i__1, i__2);
dtrevc3_((char *)"L", (char *)"B", select, n, &a[a_offset], lda, &vl[vl_offset], ldvl,
&vr[vr_offset], ldvr, n, &nout, &work[1], &c_n1, &ierr, (ftnlen)1,
(ftnlen)1);
lwork_trevc__ = (integer)work[1];
i__1 = maxwrk, i__2 = *n + lwork_trevc__;
maxwrk = max(i__1, i__2);
i__1 = maxwrk, i__2 = *n << 2;
maxwrk = max(i__1, i__2);
} else if (wantvr) {
minwrk = *n << 2;
i__1 = maxwrk,
i__2 = (*n << 1) + (*n - 1) * ilaenv_(&c__1, (char *)"DORGHR", (char *)" ", n, &c__1, n, &c_n1,
(ftnlen)6, (ftnlen)1);
maxwrk = max(i__1, i__2);
dhseqr_((char *)"S", (char *)"V", n, &c__1, n, &a[a_offset], lda, &wr[1], &wi[1], &vr[vr_offset],
ldvr, &work[1], &c_n1, info, (ftnlen)1, (ftnlen)1);
hswork = (integer)work[1];
i__1 = maxwrk, i__2 = *n + 1, i__1 = max(i__1, i__2), i__2 = *n + hswork;
maxwrk = max(i__1, i__2);
dtrevc3_((char *)"R", (char *)"B", select, n, &a[a_offset], lda, &vl[vl_offset], ldvl,
&vr[vr_offset], ldvr, n, &nout, &work[1], &c_n1, &ierr, (ftnlen)1,
(ftnlen)1);
lwork_trevc__ = (integer)work[1];
i__1 = maxwrk, i__2 = *n + lwork_trevc__;
maxwrk = max(i__1, i__2);
i__1 = maxwrk, i__2 = *n << 2;
maxwrk = max(i__1, i__2);
} else {
minwrk = *n * 3;
dhseqr_((char *)"E", (char *)"N", n, &c__1, n, &a[a_offset], lda, &wr[1], &wi[1], &vr[vr_offset],
ldvr, &work[1], &c_n1, info, (ftnlen)1, (ftnlen)1);
hswork = (integer)work[1];
i__1 = maxwrk, i__2 = *n + 1, i__1 = max(i__1, i__2), i__2 = *n + hswork;
maxwrk = max(i__1, i__2);
}
maxwrk = max(maxwrk, minwrk);
}
work[1] = (doublereal)maxwrk;
if (*lwork < minwrk && !lquery) {
*info = -13;
}
}
if (*info != 0) {
i__1 = -(*info);
xerbla_((char *)"DGEEV ", &i__1, (ftnlen)6);
return 0;
} else if (lquery) {
return 0;
}
if (*n == 0) {
return 0;
}
eps = dlamch_((char *)"P", (ftnlen)1);
smlnum = dlamch_((char *)"S", (ftnlen)1);
bignum = 1. / smlnum;
dlabad_(&smlnum, &bignum);
smlnum = sqrt(smlnum) / eps;
bignum = 1. / smlnum;
anrm = dlange_((char *)"M", n, n, &a[a_offset], lda, dum, (ftnlen)1);
scalea = FALSE_;
if (anrm > 0. && anrm < smlnum) {
scalea = TRUE_;
cscale = smlnum;
} else if (anrm > bignum) {
scalea = TRUE_;
cscale = bignum;
}
if (scalea) {
dlascl_((char *)"G", &c__0, &c__0, &anrm, &cscale, n, n, &a[a_offset], lda, &ierr, (ftnlen)1);
}
ibal = 1;
dgebal_((char *)"B", n, &a[a_offset], lda, &ilo, &ihi, &work[ibal], &ierr, (ftnlen)1);
itau = ibal + *n;
iwrk = itau + *n;
i__1 = *lwork - iwrk + 1;
dgehrd_(n, &ilo, &ihi, &a[a_offset], lda, &work[itau], &work[iwrk], &i__1, &ierr);
if (wantvl) {
*(unsigned char *)side = 'L';
dlacpy_((char *)"L", n, n, &a[a_offset], lda, &vl[vl_offset], ldvl, (ftnlen)1);
i__1 = *lwork - iwrk + 1;
dorghr_(n, &ilo, &ihi, &vl[vl_offset], ldvl, &work[itau], &work[iwrk], &i__1, &ierr);
iwrk = itau;
i__1 = *lwork - iwrk + 1;
dhseqr_((char *)"S", (char *)"V", n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vl[vl_offset], ldvl,
&work[iwrk], &i__1, info, (ftnlen)1, (ftnlen)1);
if (wantvr) {
*(unsigned char *)side = 'B';
dlacpy_((char *)"F", n, n, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, (ftnlen)1);
}
} else if (wantvr) {
*(unsigned char *)side = 'R';
dlacpy_((char *)"L", n, n, &a[a_offset], lda, &vr[vr_offset], ldvr, (ftnlen)1);
i__1 = *lwork - iwrk + 1;
dorghr_(n, &ilo, &ihi, &vr[vr_offset], ldvr, &work[itau], &work[iwrk], &i__1, &ierr);
iwrk = itau;
i__1 = *lwork - iwrk + 1;
dhseqr_((char *)"S", (char *)"V", n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vr[vr_offset], ldvr,
&work[iwrk], &i__1, info, (ftnlen)1, (ftnlen)1);
} else {
iwrk = itau;
i__1 = *lwork - iwrk + 1;
dhseqr_((char *)"E", (char *)"N", n, &ilo, &ihi, &a[a_offset], lda, &wr[1], &wi[1], &vr[vr_offset], ldvr,
&work[iwrk], &i__1, info, (ftnlen)1, (ftnlen)1);
}
if (*info != 0) {
goto L50;
}
if (wantvl || wantvr) {
i__1 = *lwork - iwrk + 1;
dtrevc3_(side, (char *)"B", select, n, &a[a_offset], lda, &vl[vl_offset], ldvl, &vr[vr_offset],
ldvr, n, &nout, &work[iwrk], &i__1, &ierr, (ftnlen)1, (ftnlen)1);
}
if (wantvl) {
dgebak_((char *)"B", (char *)"L", n, &ilo, &ihi, &work[ibal], n, &vl[vl_offset], ldvl, &ierr, (ftnlen)1,
(ftnlen)1);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
if (wi[i__] == 0.) {
scl = 1. / dnrm2_(n, &vl[i__ * vl_dim1 + 1], &c__1);
dscal_(n, &scl, &vl[i__ * vl_dim1 + 1], &c__1);
} else if (wi[i__] > 0.) {
d__1 = dnrm2_(n, &vl[i__ * vl_dim1 + 1], &c__1);
d__2 = dnrm2_(n, &vl[(i__ + 1) * vl_dim1 + 1], &c__1);
scl = 1. / dlapy2_(&d__1, &d__2);
dscal_(n, &scl, &vl[i__ * vl_dim1 + 1], &c__1);
dscal_(n, &scl, &vl[(i__ + 1) * vl_dim1 + 1], &c__1);
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
d__1 = vl[k + i__ * vl_dim1];
d__2 = vl[k + (i__ + 1) * vl_dim1];
work[iwrk + k - 1] = d__1 * d__1 + d__2 * d__2;
}
k = idamax_(n, &work[iwrk], &c__1);
dlartg_(&vl[k + i__ * vl_dim1], &vl[k + (i__ + 1) * vl_dim1], &cs, &sn, &r__);
drot_(n, &vl[i__ * vl_dim1 + 1], &c__1, &vl[(i__ + 1) * vl_dim1 + 1], &c__1, &cs,
&sn);
vl[k + (i__ + 1) * vl_dim1] = 0.;
}
}
}
if (wantvr) {
dgebak_((char *)"B", (char *)"R", n, &ilo, &ihi, &work[ibal], n, &vr[vr_offset], ldvr, &ierr, (ftnlen)1,
(ftnlen)1);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
if (wi[i__] == 0.) {
scl = 1. / dnrm2_(n, &vr[i__ * vr_dim1 + 1], &c__1);
dscal_(n, &scl, &vr[i__ * vr_dim1 + 1], &c__1);
} else if (wi[i__] > 0.) {
d__1 = dnrm2_(n, &vr[i__ * vr_dim1 + 1], &c__1);
d__2 = dnrm2_(n, &vr[(i__ + 1) * vr_dim1 + 1], &c__1);
scl = 1. / dlapy2_(&d__1, &d__2);
dscal_(n, &scl, &vr[i__ * vr_dim1 + 1], &c__1);
dscal_(n, &scl, &vr[(i__ + 1) * vr_dim1 + 1], &c__1);
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
d__1 = vr[k + i__ * vr_dim1];
d__2 = vr[k + (i__ + 1) * vr_dim1];
work[iwrk + k - 1] = d__1 * d__1 + d__2 * d__2;
}
k = idamax_(n, &work[iwrk], &c__1);
dlartg_(&vr[k + i__ * vr_dim1], &vr[k + (i__ + 1) * vr_dim1], &cs, &sn, &r__);
drot_(n, &vr[i__ * vr_dim1 + 1], &c__1, &vr[(i__ + 1) * vr_dim1 + 1], &c__1, &cs,
&sn);
vr[k + (i__ + 1) * vr_dim1] = 0.;
}
}
}
L50:
if (scalea) {
i__1 = *n - *info;
i__3 = *n - *info;
i__2 = max(i__3, 1);
dlascl_((char *)"G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wr[*info + 1], &i__2, &ierr,
(ftnlen)1);
i__1 = *n - *info;
i__3 = *n - *info;
i__2 = max(i__3, 1);
dlascl_((char *)"G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[*info + 1], &i__2, &ierr,
(ftnlen)1);
if (*info > 0) {
i__1 = ilo - 1;
dlascl_((char *)"G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wr[1], n, &ierr, (ftnlen)1);
i__1 = ilo - 1;
dlascl_((char *)"G", &c__0, &c__0, &cscale, &anrm, &i__1, &c__1, &wi[1], n, &ierr, (ftnlen)1);
}
}
work[1] = (doublereal)maxwrk;
return 0;
}
#ifdef __cplusplus
}
#endif

Some files were not shown because too many files have changed in this diff Show More