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

Merge branch 'feature-create-atoms-exclude' of github.com:erozic/lammps into create-overlap

Browse Source
This commit is contained in:
Steve Plimpton
2022-04-26 13:22:10 -06:00
parent 23d39c7d02 72011bf325
commit 782add6943
1511 changed files with 11505 additions and 40382 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
.github/CONTRIBUTING.md vendored
View File

@ -5,8 +5,8 @@ Thank your for considering to contribute to the LAMMPS software project.
The following is a set of guidelines as well as explanations of policies and work flows for contributing to the LAMMPS molecular dynamics software project. These guidelines focus on submitting issues or pull requests on the LAMMPS GitHub project.
Thus please also have a look at:
* [The guide for submitting new features in the LAMMPS manual](https://lammps.sandia.gov/doc/Modify_contribute.html)
* [The guide on programming style and requirement in the LAMMPS manual](https://lammps.sandia.gov/doc/Modify_contribute.html)
* [The guide for submitting new features in the LAMMPS manual](https://www.lammps.org/doc/Modify_contribute.html)
* [The guide on programming style and requirement in the LAMMPS manual](https://www.lammps.org/doc/Modify_style.html)
* [The GitHub tutorial in the LAMMPS manual](http://lammps.sandia.gov/doc/Howto_github.html)
## Table of Contents

23
cmake/CMakeLists.txt
View File

@ -105,8 +105,28 @@ if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
endif()
endif()
# silence excessive warnings for new Intel Compilers
if(CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM")
set(CMAKE_TUNE_DEFAULT "-Wno-tautological-constant-compare")
endif()
# silence excessive warnings for PGI/NVHPC compilers
if((CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC") OR (CMAKE_CXX_COMPILER_ID STREQUAL "PGI"))
set(CMAKE_TUNE_DEFAULT "-Minform=severe")
endif()
# silence nvcc warnings
if((PKG_KOKKOS) AND (Kokkos_ENABLE_CUDA))
set(CMAKE_TUNE_DEFAULT "${CMAKE_TUNE_DEFAULT} -Xcudafe --diag_suppress=unrecognized_pragma")
endif()
# we require C++11 without extensions. Kokkos requires at least C++14 (currently)
set(CMAKE_CXX_STANDARD 11)
if(NOT CMAKE_CXX_STANDARD)
set(CMAKE_CXX_STANDARD 11)
endif()
if(CMAKE_CXX_STANDARD LESS 11)
message(FATAL_ERROR "C++ standard must be set to at least 11")
endif()
if(PKG_KOKKOS AND (CMAKE_CXX_STANDARD LESS 14))
set(CMAKE_CXX_STANDARD 14)
endif()
@ -468,6 +488,7 @@ set(CMAKE_TUNE_FLAGS "${CMAKE_TUNE_DEFAULT}" CACHE STRING "Compiler and machine
separate_arguments(CMAKE_TUNE_FLAGS)
foreach(_FLAG ${CMAKE_TUNE_FLAGS})
target_compile_options(lammps PRIVATE ${_FLAG})
target_compile_options(lmp PRIVATE ${_FLAG})
endforeach()
########################################################################
# Basic system tests (standard libraries, headers, functions, types) #

4
cmake/Modules/Documentation.cmake
View File

@ -7,13 +7,13 @@ if(BUILD_DOC)
# Sphinx 3.x requires at least Python 3.5
if(CMAKE_VERSION VERSION_LESS 3.12)
find_package(PythonInterp 3.5 REQUIRED)
set(VIRTUALENV ${PYTHON_EXECUTABLE} -m virtualenv -p ${PYTHON_EXECUTABLE})
set(VIRTUALENV ${PYTHON_EXECUTABLE} -m venv)
else()
find_package(Python3 REQUIRED COMPONENTS Interpreter)
if(Python3_VERSION VERSION_LESS 3.5)
message(FATAL_ERROR "Python 3.5 and up is required to build the HTML documentation")
endif()
set(VIRTUALENV ${Python3_EXECUTABLE} -m virtualenv -p ${Python3_EXECUTABLE})
set(VIRTUALENV ${Python3_EXECUTABLE} -m venv)
endif()
find_package(Doxygen 1.8.10 REQUIRED)

4
cmake/Modules/Packages/COLVARS.cmake
View File

@ -19,6 +19,10 @@ endif()
add_library(colvars STATIC ${COLVARS_SOURCES})
target_compile_definitions(colvars PRIVATE -DCOLVARS_LAMMPS)
separate_arguments(CMAKE_TUNE_FLAGS)
foreach(_FLAG ${CMAKE_TUNE_FLAGS})
target_compile_options(colvars PRIVATE ${_FLAG})
endforeach()
set_target_properties(colvars PROPERTIES OUTPUT_NAME lammps_colvars${LAMMPS_MACHINE})
target_include_directories(colvars PUBLIC ${LAMMPS_LIB_SOURCE_DIR}/colvars)
# The line below is needed to locate math_eigen_impl.h

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

@ -36,3 +36,5 @@ endif()
if((CMAKE_CXX_COMPILER_ID STREQUAL "PGI") OR (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC"))
target_compile_definitions(lammps PRIVATE -DEIGEN_DONT_VECTORIZE)
endif()
target_compile_definitions(lammps PRIVATE -DEIGEN_NO_CUDA)

4
cmake/presets/kokkos-cuda.cmake
View File

@ -5,7 +5,5 @@ set(PKG_KOKKOS ON CACHE BOOL "" FORCE)
set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "" FORCE)
set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "" FORCE)
set(Kokkos_ENABLE_CUDA ON CACHE BOOL "" FORCE)
set(Kokkos_ARCH_MAXWELL50 on CACHE BOOL "" FORCE)
set(Kokkos_ARCH_PASCAL60 ON CACHE BOOL "" FORCE)
set(BUILD_OMP ON CACHE BOOL "" FORCE)
get_filename_component(NVCC_WRAPPER_CMD ${CMAKE_CURRENT_SOURCE_DIR}/../lib/kokkos/bin/nvcc_wrapper ABSOLUTE)
set(CMAKE_CXX_COMPILER ${NVCC_WRAPPER_CMD} CACHE FILEPATH "" FORCE)

2
cmake/presets/pgi.cmake
View File

@ -1,4 +1,4 @@
# preset that will enable clang/clang++ with support for MPI and OpenMP (on Linux boxes)
# preset that will enable PGI (Nvidia) compilers with support for MPI and OpenMP (on Linux boxes)
set(CMAKE_CXX_COMPILER "pgc++" CACHE STRING "" FORCE)
set(CMAKE_C_COMPILER "pgcc" CACHE STRING "" FORCE)

13
doc/Makefile
View File

@ -14,23 +14,22 @@ ANCHORCHECK = $(VENV)/bin/rst_anchor_check
SPHINXCONFIG = $(BUILDDIR)/utils/sphinx-config
MATHJAX = $(SPHINXCONFIG)/_static/mathjax
PYTHON = $(shell which python3)
DOXYGEN = $(shell which doxygen)
VIRTUALENV = virtualenv
PYTHON = $(word 3,$(shell type python3))
DOXYGEN = $(word 3,$(shell type doxygen))
HAS_PYTHON3 = NO
HAS_DOXYGEN = NO
HAS_PDFLATEX = NO
ifeq ($(shell which python3 >/dev/null 2>&1; echo $$?), 0)
ifeq ($(shell type python3 >/dev/null 2>&1; echo $$?), 0)
HAS_PYTHON3 = YES
endif
ifeq ($(shell which doxygen >/dev/null 2>&1; echo $$?), 0)
ifeq ($(shell type doxygen >/dev/null 2>&1; echo $$?), 0)
HAS_DOXYGEN = YES
endif
ifeq ($(shell which pdflatex >/dev/null 2>&1; echo $$?), 0)
ifeq ($(shell which latexmk >/dev/null 2>&1; echo $$?), 0)
ifeq ($(shell type pdflatex >/dev/null 2>&1; echo $$?), 0)
ifeq ($(shell type latexmk >/dev/null 2>&1; echo $$?), 0)
HAS_PDFLATEX = YES
endif
endif

10
doc/src/Build_development.rst
View File

@ -486,14 +486,14 @@ The following options are available.
make fix-whitespace # correct whitespace issues in files
make check-homepage # search for files with old LAMMPS homepage URLs
make fix-homepage # correct LAMMPS homepage URLs in files
make check-errordocs # search for deprecated error docs in header files
make fix-errordocs # remove error docs in header files
make check-permissions # search for files with permissions issues
make fix-permissions # correct permissions issues in files
make check # run all check targets from above
These should help to replace all TAB characters with blanks and remove
any trailing whitespace. Also all LAMMPS homepage URL references can be
updated to the location change from Sandia to the lammps.org domain.
And the permission check can remove executable permissions from non-executable
files (like source code).
These should help to make source and documentation files conforming
to some the coding style preferences of the LAMMPS developers.
Clang-format support
--------------------

1
doc/src/Commands_compute.rst
View File

@ -33,6 +33,7 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`body/local <compute_body_local>`
* :doc:`bond <compute_bond>`
* :doc:`bond/local <compute_bond_local>`
* :doc:`born/matrix <compute_born_matrix>`
* :doc:`centro/atom <compute_centro_atom>`
* :doc:`centroid/stress/atom <compute_stress_atom>`
* :doc:`chunk/atom <compute_chunk_atom>`

2
doc/src/Commands_pair.rst
View File

@ -124,7 +124,7 @@ OPT.
* :doc:`hbond/dreiding/lj (o) <pair_hbond_dreiding>`
* :doc:`hbond/dreiding/morse (o) <pair_hbond_dreiding>`
* :doc:`hdnnp <pair_hdnnp>`
* :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>`
* :doc:`ilp/graphene/hbn (t) <pair_ilp_graphene_hbn>`
* :doc:`ilp/tmd <pair_ilp_tmd>`
* :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>`
* :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>`

6
doc/src/Developer_utils.rst
View File

@ -211,6 +211,12 @@ Convenience functions
.. doxygenfunction:: logmesg(LAMMPS *lmp, const std::string &mesg)
:project: progguide
.. doxygenfunction:: errorurl
:project: progguide
.. doxygenfunction:: missing_cmd_args
:project: progguide
.. doxygenfunction:: flush_buffers(LAMMPS *lmp)
:project: progguide

1
doc/src/Errors.rst
View File

@ -11,6 +11,7 @@ them.
:maxdepth: 1
Errors_common
Errors_details
Errors_bugs
Errors_debug
Errors_messages

27
doc/src/Errors_details.rst Normal file
View File

@ -0,0 +1,27 @@
Error and warning details
=========================
Many errors or warnings are self-explanatory and thus straightforward to
resolve. However, there are also cases, where there is no single cause
and explanation, where LAMMPS can only detect symptoms of an error but
not the exact cause, or where the explanation needs to be more detailed than
what can be fit into a message printed by the program. The following are
discussions of such cases.
.. _err0001:
Unknown identifier in data file
-------------------------------
This error happens when LAMMPS encounters a line of text in an unexpected format
while reading a data file. This is most commonly cause by inconsistent header and
section data. The header section informs LAMMPS how many entries or lines are expected in the
various sections (like Atoms, Masses, Pair Coeffs, *etc.*\ ) of the data file.
If there is a mismatch, LAMMPS will either keep reading beyond the end of a section
or stop reading before the section has ended.
Such a mismatch can happen unexpectedly when the first line of the data
is *not* a comment as required by the format. That would result in
LAMMPS expecting, for instance, 0 atoms because the "atoms" header line
is treated as a comment.

39
doc/src/Howto_bpm.rst
View File

@ -3,7 +3,6 @@ Bonded particle models
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.
@ -47,33 +46,29 @@ this, LAMMPS requires :doc:`newton <newton>` bond off such that all
processors containing an atom know when a bond breaks. Additionally,
one must do either (A) or (B).
(A)
A) Use the following special bond settings
Use the following special bond settings
.. code-block:: LAMMPS
.. code-block:: LAMMPS
special_bonds lj 0 1 1 coul 1 1 1
special_bonds lj 0 1 1 coul 1 1 1
These settings accomplish two goals. First, they turn off 1-3 and 1-4
special bond lists, which are not currently supported for BPMs. As
BPMs often have dense bond networks, generating 1-3 and 1-4 special
bond lists is expensive. By setting the lj weight for 1-2 bonds to
zero, this turns off pairwise interactions. Even though there are no
charges in BPM models, setting a nonzero coul weight for 1-2 bonds
ensures all bonded neighbors are still included in the neighbor list
in case bonds break between neighbor list builds.
These settings accomplish two goals. First, they turn off 1-3 and 1-4
special bond lists, which are not currently supported for BPMs. As
BPMs often have dense bond networks, generating 1-3 and 1-4 special
bond lists is expensive. By setting the lj weight for 1-2 bonds to
zero, this turns off pairwise interactions. Even though there are no
charges in BPM models, setting a nonzero coul weight for 1-2 bonds
ensures all bonded neighbors are still included in the neighbor list
in case bonds break between neighbor list builds.
B) Alternatively, one can simply overlay pair interactions such that all
bonded particles also feel pair interactions. This can be
accomplished by using the *overlay/pair* keyword present in all bpm
bond styles and by using the following special bond settings
(B)
.. code-block:: LAMMPS
Alternatively, one can simply overlay pair interactions such that all
bonded particles also feel pair interactions. This can be accomplished
by using the *overlay/pair* keyword present in all bpm bond styles and
by using the following special bond settings
.. code-block:: LAMMPS
special_bonds lj/coul 1 1 1
special_bonds lj/coul 1 1 1
See the :doc:`Howto <Howto_broken_bonds>` page on broken bonds for
more information.

39
doc/src/Howto_elastic.rst
View File

@ -18,23 +18,52 @@ At zero temperature, it is easy to estimate these derivatives by
deforming the simulation box in one of the six directions using the
:doc:`change_box <change_box>` command and measuring the change in the
stress tensor. A general-purpose script that does this is given in the
examples/elastic directory described on the :doc:`Examples <Examples>`
examples/ELASTIC directory described on the :doc:`Examples <Examples>`
doc page.
Calculating elastic constants at finite temperature is more
challenging, because it is necessary to run a simulation that performs
time averages of differential properties. One way to do this is to
measure the change in average stress tensor in an NVT simulations when
time averages of differential properties. There are at least
3 ways to do this in LAMMPS. The most reliable way to do this is
by exploiting the relationship between elastic constants, stress
fluctuations, and the Born matrix, the second derivatives of energy
w.r.t. strain :ref:`(Ray) <Ray>`.
The Born matrix calculation has been enabled by
the :doc:`compute born/matrix <compute_born_matrix>` command,
which works for any bonded or non-bonded potential in LAMMPS.
The most expensive part of the calculation is the sampling of
the stress fluctuations. Several examples of this method are
provided in the examples/ELASTIC_T/BORN_MATRIX directory
described on the :doc:`Examples <Examples>` doc page.
A second way is to measure
the change in average stress tensor in an NVT simulations when
the cell volume undergoes a finite deformation. In order to balance
the systematic and statistical errors in this method, the magnitude of
the deformation must be chosen judiciously, and care must be taken to
fully equilibrate the deformed cell before sampling the stress
tensor. Another approach is to sample the triclinic cell fluctuations
tensor. An example of this method is provided in the
examples/ELASTIC_T/DEFORMATION directory
described on the :doc:`Examples <Examples>` doc page.
Another approach is to sample the triclinic cell fluctuations
that occur in an NPT simulation. This method can also be slow to
converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`
converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`.
We do not provide an example of this method.
A nice review of the advantages and disadvantages of all of these methods
is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier>`.
----------
.. _Ray:
**(Ray)** J. R. Ray and A. Rahman, J Chem Phys, 80, 4423 (1984).
.. _Shinoda1:
**(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
.. _Clavier:
**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).

81
doc/src/Modify_atom.rst
View File

@ -9,34 +9,34 @@ A new atom style can be created if one of the existing atom styles
does not define all the attributes you need to store and communicate
with atoms.
Atom_vec_atomic.cpp is the simplest example of an atom style.
The file ``atom_vec_atomic.cpp`` is the simplest example of an atom style.
Examining the code for others will make these instructions more clear.
Note that the :doc:`atom style hybrid <atom_style>` command can be
used to define atoms or particles which have the union of properties
of individual styles. Also the :doc:`fix property/atom <fix_property_atom>`
command can be used to add a single property (e.g. charge
or a molecule ID) to a style that does not have it. It can also be
used to add custom properties to an atom, with options to communicate
them with ghost atoms or read them from a data file. Other LAMMPS
commands can access these custom properties, as can new pair, fix,
compute styles that are written to work with these properties. For
Note that the :doc:`atom style hybrid <atom_style>` command can be used
to define atoms or particles which have the union of properties of
individual styles. Also the :doc:`fix property/atom
<fix_property_atom>` command can be used to add a single property
(e.g. charge or a molecule ID) to a style that does not have it. It can
also be used to add custom properties to an atom, with options to
communicate them with ghost atoms or read them from a data file. Other
LAMMPS commands can access these custom properties, as can new pair,
fix, compute styles that are written to work with these properties. For
example, the :doc:`set <set>` command can be used to set the values of
custom per-atom properties from an input script. All of these methods
are less work than writing code for a new atom style.
are less work than writing and testing(!) code for a new atom style.
If you follow these directions your new style will automatically work
in tandem with others via the :doc:`atom_style hybrid <atom_style>`
command.
The first step is to define a set of strings in the constructor of the
new derived class. Each string will have zero or more space-separated
variable names which are identical to those used in the atom.h header
file for per-atom properties. Note that some represent per-atom
The first step is to define a set of string lists in the constructor of
the new derived class. Each list will have zero or more comma-separated
strings that correspond to the variable names used in the ``atom.h``
header file for per-atom properties. Note that some represent per-atom
vectors (q, molecule) while other are per-atom arrays (x,v). For all
but the last 2 strings you do not need to specify any of
but the last two lists you do not need to specify any of
(id,type,x,v,f). Those are included automatically as needed in the
other strings.
other lists.
.. list-table::
@ -65,16 +65,16 @@ other strings.
* - fields_data_vel
- list of properties (in order) in the Velocities section of a data file, as read by :doc:`read_data <read_data>`
In these strings you can list variable names which LAMMPS already
defines (in some other atom style), or you can create new variable
names. You should not re-use a LAMMPS variable for something with
different meaning in your atom style. If the meaning is related, but
interpreted differently by your atom style, then using the same
variable name means a user should not use your style and the other
style together in a :doc:`atom_style hybrid <atom_style>` command.
Because there will only be one value of the variable and different
parts of LAMMPS will then likely use it differently. LAMMPS has
no way of checking for this.
In these lists you can list variable names which LAMMPS already defines
(in some other atom style), or you can create new variable names. You
should not re-use a LAMMPS variable in your atom style that is used for
something with a different meaning in another atom style. If the
meaning is related, but interpreted differently by your atom style, then
using the same variable name means a user must not use your style and
the other style together in a :doc:`atom_style hybrid <atom_style>`
command. Because there will only be one value of the variable and
different parts of LAMMPS will then likely use it differently. LAMMPS
has no way of checking for this.
If you are defining new variable names then make them descriptive and
unique to your new atom style. For example choosing "e" for energy is
@ -85,32 +85,31 @@ If any of the variable names in your new atom style do not exist in
LAMMPS, you need to add them to the src/atom.h and atom.cpp files.
Search for the word "customize" or "customization" in these 2 files to
see where to add your variable. Adding a flag to the 2nd
customization section in atom.h is only necessary if your code (e.g. a
pair style) needs to check that a per-atom property is defined. These
flags should also be set in the constructor of the atom style child
class.
see where to add your variable. Adding a flag to the 2nd customization
section in ``atom.h`` is only necessary if your code (e.g. a pair style)
needs to check that a per-atom property is defined. These flags should
also be set in the constructor of the atom style child class.
In atom.cpp, aside from the constructor and destructor, there are 3
In ``atom.cpp``, aside from the constructor and destructor, there are 3
methods that a new variable name or flag needs to be added to.
In Atom::peratom_create() when using the add_peratom() method, a
final length argument of 0 is for per-atom vectors, a length > 1 is
for per-atom arrays. Note the use of an extra per-thread flag and the
add_peratom_vary() method when last dimension of the array is
In ``Atom::peratom_create()`` when using the ``Atom::add_peratom()``
method, a cols argument of 0 is for per-atom vectors, a length >
1 is for per-atom arrays. Note the use of the extra per-thread flag and
the add_peratom_vary() method when last dimension of the array is
variable-length.
Adding the variable name to Atom::extract() enable the per-atom data
Adding the variable name to Atom::extract() enables the per-atom data
to be accessed through the :doc:`LAMMPS library interface
<Howto_library>` by a calling code, including from :doc:`Python
<Python_head>`.
The constructor of the new atom style will also typically set a few
flags which are defined at the top of atom_vec.h. If these are
flags which are defined at the top of ``atom_vec.h``. If these are
unclear, see how other atom styles use them.
The grow_pointers() method is also required to make
a copy of peratom data pointers, as explained in the code.
The grow_pointers() method is also required to make a copy of peratom
data pointers, as explained in the code.
There are a number of other optional methods which your atom style can
implement. These are only needed if you need to do something

54
doc/src/Modify_style.rst
View File

@ -223,6 +223,13 @@ and readable by all and no executable permissions. Executable
permissions (0755) should only be on shell scripts or python or similar
scripts for interpreted script languages.
You can check for these issues with the python scripts in the
:ref:`"tools/coding_standard" <coding_standard>` folder. When run
normally with a source file or a source folder as argument, they will
list all non-conforming lines. By adding the `-f` flag to the command
line, they will modify the flagged files to try removing the detected
issues.
Indentation and Placement of Braces (strongly preferred)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@ -240,6 +247,53 @@ reformatting from clang-format yields undesirable output may be
protected with placing a pair `// clang-format off` and `// clang-format
on` comments around that block.
Error or warning messages and explanations (preferred)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. versionchanged:: 27Apr2022
Starting with LAMMPS version 27 April 2022 the LAMMPS developers have
agreed on a new policy for error and warning messages.
Previously, all error and warning strings were supposed to be listed in
the class header files with an explanation. Those would then be
regularly "harvested" and transferred to alphabetically sorted lists in
the manual. To avoid excessively long lists and to reduce effort, this
came with a requirement to have rather generic error messages (e.g.
"Illegal ... command"). To identify the specific cause, the name of the
source file and the line number of the error location would be printed,
so that one could look up the cause by reading the source code.
The new policy encourages more specific error messages that ideally
indicate the cause directly and no further lookup would be needed.
This is aided by using the `{fmt} library <https://fmt.dev>`_ to convert
the Error class commands so that they take a variable number of arguments
and error text will be treated like a {fmt} syntax format string.
Error messages should still kept to a single line or two lines at the most.
For more complex explanations or errors that have multiple possible
reasons, a paragraph should be added to the `Error_details` page with an
error code reference (e.g. ``.. _err0001:``) then the utility function
:cpp:func:`utils::errorurl() <LAMMPS_NS::utils::errorurl>` can be used
to generate an URL that will directly lead to that paragraph. An error
for missing arguments can be easily generated using the
:cpp:func:`utils::missing_cmd_args()
<LAMMPS_NS::utils::missing_cmd_args>` convenience function.
The transformation of existing LAMMPS code to this new scheme is ongoing
and - given the size of the LAMMPS source code - will take a significant
amount of time until completion. However, for new code following the
new approach is strongly preferred. The expectation is that the new
scheme will make it easier for LAMMPS users, developers, and
maintainers.
An example for this approach would be the
``src/read_data.cpp`` and ``src/atom.cpp`` files that implement the
:doc:`read_data <read_data>` and :doc:`atom_modify <atom_modify>`
commands and that may create :ref:`"Unknown identifier in data file" <err0001>`
errors that seem difficult to debug for users because they may have
one of multiple possible reasons, and thus require some additional explanations.
Programming language standards (required)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

35
doc/src/Tools.rst
View File

@ -87,7 +87,7 @@ Miscellaneous tools
.. table_from_list::
:columns: 6
* :ref:`CMake <cmake>`
* :ref:`LAMMPS coding standards <coding_standard>`
* :ref:`emacs <emacs>`
* :ref:`i-pi <ipi>`
* :ref:`kate <kate>`
@ -189,27 +189,32 @@ for the :doc:`chain benchmark <Speed_bench>`.
----------
.. _cmake:
.. _coding_standard:
CMake tools
-----------
LAMMPS coding standard
----------------------
The ``cmbuild`` script is a wrapper around using ``cmake --build <dir>
--target`` and allows compiling LAMMPS in a :ref:`CMake build folder
<cmake_build>` with a make-like syntax regardless of the actual build
tool and the specific name of the program used (e.g. ``ninja-v1.10`` or
``gmake``) when using ``-D CMAKE_MAKE_PROGRAM=<name>``.
The ``coding_standard`` folder contains multiple python scripts to
check for and apply some LAMMPS coding conventions. The following
scripts are available:
.. parsed-literal::
Usage: cmbuild [-v] [-h] [-C <dir>] [-j <num>] [<target>]
permissions.py # detects if sources have executable permissions and scripts have not
whitespace.py # detects TAB characters and trailing whitespace
homepage.py # detects outdated LAMMPS homepage URLs (pointing to sandia.gov instead of lammps.org)
errordocs.py # detects deprecated error docs in header files
Options:
-h print this message
-j <NUM> allow processing of NUM concurrent tasks
-C DIRECTORY execute build in folder DIRECTORY
-v produce verbose output
The tools need to be given the main folder of the LAMMPS distribution
or individual file names as argument and will by default check them
and report any non-compliance. With the optional ``-f`` argument the
corresponding script will try to change the non-compliant file(s) to
match the conventions.
For convenience this scripts can also be invoked by the make file in
the ``src`` folder with, `make check-whitespace` or `make fix-whitespace`
to either detect or edit the files. Correspondingly for the other python
scripts. `make check` will run all checks.
----------

1
doc/src/compute.rst
View File

@ -179,6 +179,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
* :doc:`body/local <compute_body_local>` - attributes of body sub-particles
* :doc:`bond <compute_bond>` - energy of each bond sub-style
* :doc:`bond/local <compute_bond_local>` - distance and energy of each bond
* :doc:`born/matrix <compute_born_matrix>` - second derivative or potential with respect to strain
* :doc:`centro/atom <compute_centro_atom>` - centro-symmetry parameter for each atom
* :doc:`centroid/stress/atom <compute_stress_atom>` - centroid based stress tensor for each atom
* :doc:`chunk/atom <compute_chunk_atom>` - assign chunk IDs to each atom

213
doc/src/compute_born_matrix.rst Normal file
View File

@ -0,0 +1,213 @@
.. index:: compute born/matrix
compute born/matrix command
===========================
Syntax
""""""
.. parsed-literal::
compute ID group-ID born/matrix keyword value ...
* ID, group-ID are documented in :doc:`compute <compute>` command
* born/matrix = style name of this compute command
* zero or more keyword/value pairs may be appended
.. parsed-literal::
keyword = *numdiff*
*numdiff* values = delta virial-ID
delta = magnitude of strain (dimensionless)
virial-ID = ID of pressure compute for virial (string)
Examples
""""""""
.. code-block:: LAMMPS
compute 1 all born/matrix
compute 1 all born/matrix bond angle
compute 1 all born/matrix numdiff 1.0e-4 myvirial
Description
"""""""""""
Define a compute that calculates
:math:`\frac{\partial{}^2U}{\partial\varepsilon_{i}\partial\varepsilon_{j}}` the
second derivatives of the potential energy :math:`U` w.r.t. strain
tensor :math:`\varepsilon` elements. These values are related to:
.. math::
C^{B}_{i,j}=\frac{1}{V}\frac{\partial{}^2U}{\partial{}\varepsilon_{i}\partial\varepsilon_{j}}
also called the Born term of elastic constants in the stress-stress fluctuation
formalism. This quantity can be used to compute the elastic constant tensor.
Using the symmetric Voigt notation, the elastic constant tensor can be written
as a 6x6 symmetric matrix:
.. math::
C_{i,j} = \langle{}C^{B}_{i,j}\rangle
+ \frac{V}{k_{B}T}\left(\langle\sigma_{i}\sigma_{j}\rangle\right.
\left.- \langle\sigma_{i}\rangle\langle\sigma_{j}\rangle\right)
+ \frac{Nk_{B}T}{V}
\left(\delta_{i,j}+(\delta_{1,i}+\delta_{2,i}+\delta_{3,i})\right.
\left.*(\delta_{1,j}+\delta_{2,j}+\delta_{3,j})\right)
In the above expression, :math:`\sigma` stands for the virial stress
tensor, :math:`\delta` is the Kronecker delta and the usual notation apply for
the number of particle, the temperature and volume respectively :math:`N`,
:math:`T` and :math:`V`. :math:`k_{B}` is the Boltzmann constant.
The Born term is a symmetric 6x6 matrix, as is the matrix of second derivatives
of potential energy w.r.t strain,
whose 21 independent elements are output in this order:
.. math::
\begin{matrix}
C_{1} & C_{7} & C_{8} & C_{9} & C_{10} & C_{11} \\
C_{7} & C_{2} & C_{12} & C_{13} & C_{14} & C_{15} \\
\vdots & C_{12} & C_{3} & C_{16} & C_{17} & C_{18} \\
\vdots & C_{13} & C_{16} & C_{4} & C_{19} & C_{20} \\
\vdots & \vdots & \vdots & C_{19} & C_{5} & C_{21} \\
\vdots & \vdots & \vdots & \vdots & C_{21} & C_{6}
\end{matrix}
in this matrix the indices of :math:`C_{k}` value are the corresponding element
:math:`k` in the global vector output by this compute. Each term comes from the sum
of the derivatives of every contribution to the potential energy
in the system as explained in :ref:`(VanWorkum)
<VanWorkum>`.
The output can be accessed using usual Lammps routines:
.. code-block:: LAMMPS
compute 1 all born/matrix
compute 2 all pressure NULL virial
variable S1 equal -c_2[1]
variable S2 equal -c_2[2]
variable S3 equal -c_2[3]
variable S4 equal -c_2[4]
variable S5 equal -c_2[5]
variable S6 equal -c_2[6]
fix 1 all ave/time 1 1 1 v_S1 v_S2 v_S3 v_S4 v_S5 v_S6 c_1[*] file born.out
In this example, the file *born.out* will contain the information needed to
compute the first and second terms of the elastic constant matrix in a post
processing procedure. The other required quantities can be accessed using any
other *LAMMPS* usual method. Several examples of this method are
provided in the examples/ELASTIC_T/BORN_MATRIX directory
described on the :doc:`Examples <Examples>` doc page.
NOTE: In the above :math:`C_{i,j}` computation, the fluctuation
term involving the virial stress tensor :math:`\sigma` is the
covariance between each elements. In a
solid the stress fluctuations can vary rapidly, while average
fluctuations can be slow to converge.
A detailed analysis of the convergence rate of all the terms in
the elastic tensor
is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier2>`.
Two different computation methods for the Born matrix are implemented in this
compute and are mutually exclusive.
The first one is a direct computation from the analytical formula from the
different terms of the potential used for the simulations :ref:`(VanWorkum)
<VanWorkum>`. However, the implementation of such derivations must be done
for every potential form. This has not been done yet and can be very
complicated for complex potentials. At the moment a warning message is
displayed for every term that is not supporting the compute at the moment.
This method is the default for now.
The second method uses finite differences of energy to numerically approximate
the second derivatives :ref:`(Zhen) <Zhen>`. This is useful when using
interaction styles for which the analytical second derivatives have not been
implemented. In this cases, the compute applies linear strain fields of
magnitude *delta* to all the atoms relative to a point at the center of the
box. The strain fields are in six different directions, corresponding to the
six Cartesian components of the stress tensor defined by LAMMPS. For each
direction it applies the strain field in both the positive and negative senses,
and the new stress virial tensor of the entire system is calculated after each.
The difference in these two virials divided by two times *delta*, approximates
the corresponding components of the second derivative, after applying a
suitable unit conversion.
.. note::
It is important to choose a suitable value for delta, the magnitude of
strains that are used to generate finite difference
approximations to the exact virial stress. For typical systems, a value in
the range of 1 part in 1e5 to 1e6 will be sufficient.
However, the best value will depend on a multitude of factors
including the stiffness of the interatomic potential, the thermodynamic
state of the material being probed, and so on. The only way to be sure
that you have made a good choice is to do a sensitivity study on a
representative atomic configuration, sweeping over a wide range of
values of delta. If delta is too small, the output values will vary
erratically due to truncation effects. If delta is increased beyond a
certain point, the output values will start to vary smoothly with
delta, due to growing contributions from higher order derivatives. In
between these two limits, the numerical virial values should be largely
independent of delta.
The keyword requires the additional arguments *delta* and *virial-ID*.
*delta* gives the size of the applied strains. *virial-ID* gives
the ID string of the pressure compute that provides the virial stress tensor,
requiring that it use the virial keyword e.g.
.. code-block:: LAMMPS
compute myvirial all pressure NULL virial
compute 1 all born/matrix numdiff 1.0e-4 myvirial
**Output info:**
This compute calculates a global vector with 21 values that are
the second derivatives of the potential energy w.r.t. strain.
The values are in energy units.
The values are ordered as explained above. These values can be used
by any command that uses global values from a compute as input. See
the :doc:`Howto output <Howto_output>` doc page for an overview of
LAMMPS output options.
The array values calculated by this compute are all "extensive".
Restrictions
""""""""""""
This compute is part of the EXTRA-COMPUTE package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info. LAMMPS was built with that package. See
the :doc:`Build package <Build_package>` page for more info.
The Born term can be decomposed as a product of two terms. The first one is a
general term which depends on the configuration. The second one is specific to
every interaction composing your force field (non-bonded, bonds, angle...).
Currently not all LAMMPS interaction styles implement the *born_matrix* method
giving first and second order derivatives and LAMMPS will exit with an error if
this compute is used with such interactions unless the *numdiff* option is
also used. The *numdiff* option cannot be used with any other keyword. In this
situation, LAMMPS will also exit with an error.
Default
"""""""
none
----------
.. _VanWorkum:
**(Van Workum)** K. Van Workum et al., J. Chem. Phys. 125 144506 (2006)
.. _Clavier2:
**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).
.. _Zhen:
**(Zhen)** Y. Zhen, C. Chu, Computer Physics Communications 183(2012)261-265

7
doc/src/compute_rigid_local.rst
View File

@ -127,19 +127,16 @@ The *vx*, *vy*, *vz*, *fx*, *fy*, *fz* attributes are components of
the COM velocity and force on the COM of the body.
The *omegax*, *omegay*, and *omegaz* attributes are the angular
velocity components of the body around its COM.
velocity components of the body in the system frame around its COM.
The *angmomx*, *angmomy*, and *angmomz* attributes are the angular
momentum components of the body around its COM.
momentum components of the body in the system frame around its COM.
The *quatw*, *quati*, *quatj*, and *quatk* attributes are the
components of the 4-vector quaternion representing the orientation of
the rigid body. See the :doc:`set <set>` command for an explanation of
the quaternion vector.
The *angmomx*, *angmomy*, and *angmomz* attributes are the angular
momentum components of the body around its COM.
The *tqx*, *tqy*, *tqz* attributes are components of the torque acting
on the body around its COM.

45
doc/src/compute_temp_profile.rst
View File

@ -76,21 +76,28 @@ velocity for each atom. Note that if there is only one atom in the
bin, its thermal velocity will thus be 0.0.
After the spatially-averaged velocity field has been subtracted from
each atom, the temperature is calculated by the formula KE = (dim\*N
- dim\*Nx\*Ny\*Nz) k T/2, where KE = total kinetic energy of the group of
atoms (sum of 1/2 m v\^2), dim = 2 or 3 = dimensionality of the
simulation, N = number of atoms in the group, k = Boltzmann constant,
and T = temperature. The dim\*Nx\*Ny\*Nz term are degrees of freedom
subtracted to adjust for the removal of the center-of-mass velocity in
each of Nx\*Ny\*Nz bins, as discussed in the :ref:`(Evans) <Evans1>` paper.
each atom, the temperature is calculated by the formula
*KE* = (*dim\*N* - *Ns\*Nx\*Ny\*Nz* - *extra* ) *k* *T*/2, where *KE* = total
kinetic energy of the group of atoms (sum of 1/2 *m* *v*\^2), *dim* = 2
or 3 = dimensionality of the simulation, *Ns* = 0, 1, 2 or 3 for
streaming velocity subtracted in 0, 1, 2 or 3 dimensions, *extra* = extra
degrees-of-freedom, *N* = number of atoms in the group, *k* = Boltzmann
constant, and *T* = temperature. The *Ns\*Nx\*Ny\*Nz* term is degrees
of freedom subtracted to adjust for the removal of the center-of-mass
velocity in each direction of the *Nx\*Ny\*Nz* bins, as discussed in the
:ref:`(Evans) <Evans1>` paper. The extra term defaults to (*dim* - *Ns*)
and accounts for overall conservation of center-of-mass velocity across
the group in directions where streaming velocity is *not* subtracted. This
can be altered using the *extra* option of the
:doc:`compute_modify <compute_modify>` command.
If the *out* keyword is used with a *tensor* value, which is the
default, a kinetic energy tensor, stored as a 6-element vector, is
also calculated by this compute for use in the computation of a
pressure tensor. The formula for the components of the tensor is the
same as the above formula, except that v\^2 is replaced by vx\*vy for
the xy component, etc. The 6 components of the vector are ordered xx,
yy, zz, xy, xz, yz.
same as the above formula, except that *v*\^2 is replaced by *vx\*vy* for
the xy component, etc. The 6 components of the vector are ordered *xx,
yy, zz, xy, xz, yz.*
If the *out* keyword is used with a *bin* value, the count of atoms
and computed temperature for each bin are stored for output, as an
@ -123,10 +130,20 @@ needed, the subtracted degrees-of-freedom can be altered using the
.. note::
When using the *out* keyword with a value of *bin*, the
calculated temperature for each bin does not include the
degrees-of-freedom adjustment described in the preceding paragraph,
for fixes that constrain molecular motion. It does include the
adjustment due to the *extra* option, which is applied to each bin.
calculated temperature for each bin includes the degrees-of-freedom
adjustment described in the preceding paragraph for fixes that
constrain molecular motion, as well as the adjustment due to
the *extra* option (which defaults to *dim* - *Ns* as described above),
by fractionally applying them based on the fraction of atoms in each
bin. As a result, the bin degrees-of-freedom summed over all bins exactly
equals the degrees-of-freedom used in the scalar temperature calculation,
:math:`\Sigma N_{DOF_i} = N_{DOF}` and the corresponding relation for temperature
is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T`.
These relations will breakdown in cases where the adjustment
exceeds the actual number of degrees-of-freedom in a bin. This could happen
if a bin is empty or in situations where rigid molecules
are non-uniformly distributed, in which case the reported
temperature within a bin may not be accurate.
See the :doc:`Howto thermostat <Howto_thermostat>` page for a
discussion of different ways to compute temperature and perform

10
doc/src/create_atoms.rst
View File

@ -314,6 +314,16 @@ radius is taken into account so that all new molecules will be created
at locations not closer than (*radius* + molecule radius) from the location
of any existing atom in the system.
.. note::
Checking for overlaps is a very costly operation (O(N) for each new atom/molecule,
where N is the number of existing atoms) and the intended use of this keyword is,
for example, adding small amounts of new atoms/molecules to relatively sparse systems
mid simulation (between consecutive runs), i.e. where running an energy minimization
procedure isn't an option.
In any case, the use of the *maxtries* keyword in combination with *overlap* is
highly recommended.
The *units* keyword determines the meaning of the distance units used
to specify the coordinates of the one particle created by the *single*
style. A *box* value selects standard distance units as defined by

87
doc/src/fix_adapt.rst
View File

@ -14,7 +14,7 @@ Syntax
* adapt = style name of this fix command
* N = adapt simulation settings every this many timesteps
* one or more attribute/arg pairs may be appended
* attribute = *pair* or *bond* or *kspace* or *atom*
* attribute = *pair* or *bond* or *angle* or *kspace* or *atom*
.. parsed-literal::
@ -28,11 +28,16 @@ Syntax
bparam = parameter to adapt over time
I = type bond to set parameter for
v_name = variable with name that calculates value of bparam
*angle* args = astyle aparam I v_name
astyle = angle style name, e.g. harmonic
aparam = parameter to adapt over time
I = type angle to set parameter for
v_name = variable with name that calculates value of aparam
*kspace* arg = v_name
v_name = variable with name that calculates scale factor on K-space terms
*atom* args = aparam v_name
aparam = parameter to adapt over time
v_name = variable with name that calculates value of aparam
*atom* args = atomparam v_name
atomparam = parameter to adapt over time
v_name = variable with name that calculates value of atomparam
* zero or more keyword/value pairs may be appended
* keyword = *scale* or *reset* or *mass*
@ -283,30 +288,62 @@ operates. The only difference is that now a bond coefficient for a
given bond type is adapted.
A wild-card asterisk can be used in place of or in conjunction with
the bond type argument to set the coefficients for multiple bond types.
This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N = the number of
atom types, then an asterisk with no numeric values means all types
from 1 to N. A leading asterisk means all types from 1 to n (inclusive).
A trailing asterisk means all types from n to N (inclusive). A middle
asterisk means all types from m to n (inclusive).
the bond type argument to set the coefficients for multiple bond
types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N =
the number of bond types, then an asterisk with no numeric values
means all types from 1 to N. A leading asterisk means all types from
1 to n (inclusive). A trailing asterisk means all types from n to N
(inclusive). A middle asterisk means all types from m to n
(inclusive).
Currently *bond* does not support bond_style hybrid nor bond_style
hybrid/overlay as bond styles. The only bonds that currently are
working with fix_adapt are
hybrid/overlay as bond styles. The bond styles that currently work
with fix_adapt are
+------------------------------------+-------+------------+
| :doc:`class2 <bond_class2>` | r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`fene <bond_fene>` | k, r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`gromos <bond_gromos>` | k, r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`harmonic <bond_harmonic>` | k,r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`morse <bond_morse>` | r0 | type bonds |
+------------------------------------+-------+------------+
| :doc:`nonlinear <bond_nonlinear>` | r0 | type bonds |
+------------------------------------+-------+------------+
+------------------------------------+-------+-----------------+
| :doc:`class2 <bond_class2>` | r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`fene <bond_fene>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`fene/nm <bond_fene>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`gromos <bond_gromos>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`harmonic <bond_harmonic>` | k,r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`morse <bond_morse>` | r0 | type bonds |
+------------------------------------+-------+-----------------+
| :doc:`nonlinear <bond_nonlinear>` | epsilon,r0 | type bonds |
+------------------------------------+-------+-----------------+
----------
The *angle* keyword uses the specified variable to change the value of
an angle coefficient over time, very similar to how the *pair* keyword
operates. The only difference is that now an angle coefficient for a
given angle type is adapted.
A wild-card asterisk can be used in place of or in conjunction with
the angle type argument to set the coefficients for multiple angle
types. This takes the form "\*" or "\*n" or "n\*" or "m\*n". If N =
the number of angle types, then an asterisk with no numeric values
means all types from 1 to N. A leading asterisk means all types from
1 to n (inclusive). A trailing asterisk means all types from n to N
(inclusive). A middle asterisk means all types from m to n
(inclusive).
Currently *angle* does not support angle_style hybrid nor angle_style
hybrid/overlay as angle styles. The angle styles that currently work
with fix_adapt are
+------------------------------------+-------+-----------------+
| :doc:`harmonic <angle_harmonic>` | k,theta0 | type angles |
+------------------------------------+-------+-----------------+
| :doc:`cosine <angle_cosine>` | k | 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.
----------

6
doc/src/fix_nve.rst
View File

@ -35,6 +35,10 @@ consistent with the microcanonical ensemble (NVE) provided there
are (full) periodic boundary conditions and no other "manipulations"
of the system (e.g. fixes that modify forces or velocities).
This fix invokes the velocity form of the
Stoermer-Verlet time integration algorithm (velocity-Verlet). Other
time integration options can be invoked using the :doc:`run_style <run_style>` command.
----------
.. include:: accel_styles.rst
@ -57,7 +61,7 @@ Restrictions
Related commands
""""""""""""""""
:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`
:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`, :doc:`run_style <run_style>`
Default
"""""""

16
doc/src/fix_property_atom.rst
View File

@ -304,13 +304,15 @@ uninterrupted fashion.
.. warning::
When reading data from a restart file, this fix command has to be
specified **exactly** the same was in the input script that created
the restart file. LAMMPS will only check whether a fix is of the
same style and has the same fix ID and in case of a match will then
try to initialize the fix with the data stored in the binary
restart file. If the names and associated date types in the new
fix property/atom command do not match the old one exactly, data
can be corrupted or LAMMPS may crash.
specified **after** the *read_restart* command and **exactly** the
same was in the input script that created the restart file. LAMMPS
will only check whether a fix is of the same style and has the same
fix ID and in case of a match will then try to initialize the fix
with the data stored in the binary restart file. If the names and
associated date types in the new fix property/atom command do not
match the old one exactly, data can be corrupted or LAMMPS may crash.
If the fix is specified **before** the *read_restart* command its
data will not be restored.
None of the :doc:`fix_modify <fix_modify>` options are relevant to
this fix. No global or per-atom quantities are stored by this fix for

2
doc/src/fix_ttm.rst
View File

@ -217,7 +217,7 @@ units used.
.. note::
The electronic temperature at each grid point must be a non-zero
positive value, both initially, and as the temperature evovles over
positive value, both initially, and as the temperature evolves over
time. Thus you must use either the *set* or *infile* keyword or be
restarting a simulation that used this fix previously.

16
doc/src/group.rst
View File

@ -258,11 +258,17 @@ assignment is made at the beginning of the minimization, but not
during the iterations of the minimizer.
The point in the timestep at which atoms are assigned to a dynamic
group is after the initial stage of velocity Verlet time integration
has been performed, and before neighbor lists or forces are computed.
This is the point in the timestep where atom positions have just
changed due to the time integration, so the region criterion should be
accurate, if applied.
group is after interatomic forces have been computed, but before any
fixes which alter forces or otherwise update the system have been
invoked. This means that atom positions have been updated, neighbor
lists and ghost atoms are current, and both intermolecular and
intramolecular forces have been calculated based on the new
coordinates. Thus the region criterion, if applied, should be
accurate. Also, any computes invoked by an atom-style variable should
use updated information for that timestep, e.g. potential energy/atom
or coordination number/atom. Similarly, fixes or computes which are
invoked after that point in the timestep, should operate on the new
group of atoms.
.. note::

7
doc/src/pair_ilp_graphene_hbn.rst
View File

@ -1,8 +1,11 @@
.. index:: pair_style ilp/graphene/hbn
.. index:: pair_style ilp/graphene/hbn/opt
pair_style ilp/graphene/hbn command
===================================
Accelerator Variant: *ilp/graphene/hbn/opt*
Syntax
""""""
@ -125,6 +128,10 @@ headings) the following commands could be included in an input script:
----------
.. include:: accel_styles.rst
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

3
doc/src/run_style.rst
View File

@ -67,7 +67,8 @@ Description
Choose the style of time integrator used for molecular dynamics
simulations performed by LAMMPS.
The *verlet* style is a standard velocity-Verlet integrator.
The *verlet* style is the velocity form of the
Stoermer-Verlet time integration algorithm (velocity-Verlet)
----------

2
doc/src/thermo_modify.rst
View File

@ -252,6 +252,6 @@ flush = no, and temp/press = compute IDs defined by thermo_style.
The defaults for the line and format options depend on the thermo style.
For styles "one" and "custom", the line and format defaults are "one",
"%10d", and "%12.8g". For style "multi", the line and format defaults
"%10d", and "%14.8g". For style "multi", the line and format defaults
are "multi", "%14d", and "%14.4f". For style "yaml", the line and format
defaults are "%d" and "%.15g".

2
doc/src/thermo_style.rst
View File

@ -10,7 +10,7 @@ Syntax
thermo_style style args
* style = *one* or *multi* *yaml* or *custom*
* style = *one* or *multi* or *yaml* or *custom*
* args = list of arguments for a particular style
.. parsed-literal::

1
doc/utils/requirements.txt
View File

@ -6,3 +6,4 @@ breathe
Pygments
six
pyyaml
wheel

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

@ -315,6 +315,7 @@ borophene
Botero
Botu
Bouguet
Bourasseau
Bourne
boxcolor
boxhi
@ -342,6 +343,7 @@ Broglie
brownian
brownw
Broyden
Brusselle
Bryantsev
Btarget
btype
@ -674,6 +676,7 @@ Deresiewicz
Derjagin
Derjaguin
Derlet
Desbiens
Deserno
Destree
destructor
@ -787,6 +790,7 @@ Dullweber
dumpfile
Dunbrack
Dunweg
Dupend
Dupont
dUs
dV
@ -1669,6 +1673,7 @@ Kusters
Kutta
Kuznetsov
kx
Lachet
Lackmann
Ladd
lagrangian
@ -3188,6 +3193,7 @@ stochastically
stochasticity
Stockmayer
Stoddard
Stoermer
stoichiometric
stoichiometry
Stokesian
@ -3599,6 +3605,7 @@ Voronoi
VORONOI
Vorselaars
Voth
Voyiatzis
vpz
vratio
Vries
@ -3665,6 +3672,7 @@ wn
Wolde
workflow
workflows
Workum
Worley
Wriggers
Wuppertal

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/born_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
3.36316 1.87373 1.87607 -0.00346 -0.00172 -0.00104
1.87373 3.36170 1.87425 0.00443 0.00033 0.00014
1.87607 1.87425 3.36573 0.00143 0.00155 0.00127
-0.00346 0.00443 0.00143 1.87425 0.00127 0.00033
-0.00172 0.00033 0.00155 0.00127 1.87607 -0.00346
-0.00104 0.00014 0.00127 0.00033 -0.00346 1.87373

118
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py Normal file
View File

@ -0,0 +1,118 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import sys
import numpy as np
def reduce_Born(Cf):
C = np.zeros((6,6), dtype=np.float64)
C[0,0] = Cf[0]
C[1,1] = Cf[1]
C[2,2] = Cf[2]
C[3,3] = Cf[3]
C[4,4] = Cf[4]
C[5,5] = Cf[5]
C[0,1] = Cf[6]
C[0,2] = Cf[7]
C[0,3] = Cf[8]
C[0,4] = Cf[9]
C[0,5] = Cf[10]
C[1,2] = Cf[11]
C[1,3] = Cf[12]
C[1,4] = Cf[13]
C[1,5] = Cf[14]
C[2,3] = Cf[15]
C[2,4] = Cf[16]
C[2,5] = Cf[17]
C[3,4] = Cf[18]
C[3,5] = Cf[19]
C[4,5] = Cf[20]
C = np.where(C,C,C.T)
return C
def compute_delta():
D = np.zeros((3,3,3,3))
for a in range(3):
for b in range(3):
for m in range(3):
for n in range(3):
D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
d = np.zeros((6,6))
d[0,0] = D[0,0,0,0]
d[1,1] = D[1,1,1,1]
d[2,2] = D[2,2,2,2]
d[3,3] = D[1,2,1,2]
d[4,4] = D[0,2,0,2]
d[5,5] = D[0,1,0,1]
d[0,1] = D[0,0,1,1]
d[0,2] = D[0,0,2,2]
d[0,3] = D[0,0,1,2]
d[0,4] = D[0,0,0,2]
d[0,5] = D[0,0,0,1]
d[1,2] = D[1,1,2,2]
d[1,3] = D[1,1,1,2]
d[1,4] = D[1,1,0,2]
d[1,5] = D[1,1,0,1]
d[2,3] = D[2,2,1,2]
d[2,4] = D[2,2,0,2]
d[2,5] = D[2,2,0,1]
d[3,4] = D[1,2,0,2]
d[3,5] = D[1,2,0,1]
d[4,5] = D[0,2,0,1]
d = np.where(d,d,d.T)
return d
def write_matrix(C, filename):
with open(filename, 'w') as f:
f.write("# Cij Matrix from post process computation\n")
for i in C:
f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
)
)
return
def main():
N = 500
vol = 27.047271**3 * 10**-30 # m^3
T = 60 # K
kb = 1.380649 * 10**-23 # J/K
kbT = T*kb # J
kcalmol2J = 4183.9954/(6.022*10**23)
born = np.loadtxt('born.out')
stre = np.loadtxt('vir.out')
stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa
try:
mean_born = np.mean(born[:, 1:], axis=0)
except IndexError:
mean_born = born[1:]
CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa
Cs = vol/kbT*np.cov(stre[:,1:].T)
Ct = N*kbT/vol * compute_delta()
C = CB - Cs + Ct
write_matrix(CB, 'born_matrix.out')
write_matrix(Cs, 'stre_matrix.out')
write_matrix(Ct, 'temp_matrix.out')
write_matrix(C, 'full_matrix.out')
C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
print(C*10**-9)
print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
return
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
raise SystemExit("User interruption.")

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/full_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
2.18161 1.13726 1.16596 -0.01607 -0.02637 0.00291
1.13726 2.20242 1.16714 0.00386 -0.05820 0.02644
1.16596 1.16714 2.24704 -0.00354 -0.00368 0.02714
-0.01607 0.00386 -0.00354 1.43706 0.00210 0.01003
-0.02637 -0.05820 -0.00368 0.00210 1.37530 0.01401
0.00291 0.02644 0.02714 0.01003 0.01401 1.42403

154
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov Normal file
View File

@ -0,0 +1,154 @@
# Analytical calculation
# of Born matrix
# Note that because of cubic symmetry and central forces, we have:
# C11, pure axial == positive mean value: 1,2,3
# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
# adjustable parameters
units real
variable nsteps index 100000 # length of run
variable nthermo index 1000 # thermo output interval
variable nlat equal 5 # size of box
variable T equal 60. # Temperature in K
variable rho equal 5.405 # Lattice spacing in A
atom_style atomic
lattice fcc ${rho}
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
create_box 1 box
create_atoms 1 box
mass * 39.948
velocity all create ${T} 87287 loop geom
velocity all zero linear
pair_style lj/cut 12.0
pair_coeff 1 1 0.238067 3.405
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
variable vol equal vol
thermo 100
fix aL all ave/time 1 1 1 v_vol ave running
fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
run 20000
unfix NPT
variable newL equal "f_aL^(1./3.)"
change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
unfix aL
reset_timestep 0
# Conversion variables
variable kb equal 1.38065e-23 # J/K
variable Myvol equal "vol*10^-30" # Volume in m^3
variable kbt equal "v_kb*v_T"
variable Nat equal atoms
variable Rhokbt equal "v_kbt*v_Nat/v_Myvol"
variable at2Pa equal 101325
variable kcalmol2J equal "4183.9954/(6.022e23)"
variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms
variable Pa2GPa equal 1e-9
# Born compute giving <C^b> terms
compute born all born/matrix
# The six virial stress component to compute <C^fl>
compute VIR all pressure NULL virial
variable s1 equal "-c_VIR[1]*v_at2Pa"
variable s2 equal "-c_VIR[2]*v_at2Pa"
variable s3 equal "-c_VIR[3]*v_at2Pa"
variable s6 equal "-c_VIR[4]*v_at2Pa"
variable s5 equal "-c_VIR[5]*v_at2Pa"
variable s4 equal "-c_VIR[6]*v_at2Pa"
variable press equal press
# Average of Born term and vector to store stress
# for post processing
fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
fix APR all ave/time 1 1 1 v_press ave running
fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
thermo ${nthermo}
thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
thermo_modify line multi
fix 1 all nvt temp $T $T 100
run ${nsteps}
# Compute vector averages
# Note the indice switch.
# LAMMPS convention is NOT the Voigt notation.
variable aves1 equal "ave(f_VEC[1])"
variable aves2 equal "ave(f_VEC[2])"
variable aves3 equal "ave(f_VEC[3])"
variable aves4 equal "ave(f_VEC[6])"
variable aves5 equal "ave(f_VEC[5])"
variable aves6 equal "ave(f_VEC[4])"
# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
# is numerically instable. Here we go through the <(s-<s>)^2>
# definition.
# Computing difference relative to average values
variable ds1 vector "f_VEC[1]-v_aves1"
variable ds2 vector "f_VEC[2]-v_aves2"
variable ds3 vector "f_VEC[3]-v_aves3"
variable ds4 vector "f_VEC[4]-v_aves4"
variable ds5 vector "f_VEC[5]-v_aves5"
variable ds6 vector "f_VEC[6]-v_aves6"
# Squaring and averaging
variable dds1 vector "v_ds1*v_ds1"
variable dds2 vector "v_ds2*v_ds2"
variable dds3 vector "v_ds3*v_ds3"
variable vars1 equal "ave(v_dds1)"
variable vars2 equal "ave(v_dds2)"
variable vars3 equal "ave(v_dds3)"
variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
variable dds12 vector "v_ds1*v_ds2"
variable dds13 vector "v_ds1*v_ds3"
variable dds23 vector "v_ds2*v_ds3"
variable vars12 equal "ave(v_dds12)"
variable vars13 equal "ave(v_dds13)"
variable vars23 equal "ave(v_dds23)"
variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
variable dds4 vector "v_ds4*v_ds4"
variable dds5 vector "v_ds5*v_ds5"
variable dds6 vector "v_ds6*v_ds6"
variable vars4 equal "ave(v_dds4)"
variable vars5 equal "ave(v_dds5)"
variable vars6 equal "ave(v_dds6)"
variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
print """
C11 = ${aC11}
C12 = ${aC12}
C44 = ${aC44}
"""

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/stre_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
1.22342 0.73647 0.71011 0.01261 0.02465 -0.00395
0.73647 1.20115 0.70711 0.00057 0.05854 -0.02630
0.71011 0.70711 1.16055 0.00497 0.00524 -0.02587
0.01261 0.00057 0.00497 0.45813 -0.00083 -0.00969
0.02465 0.05854 0.00524 -0.00083 0.52170 -0.01747
-0.00395 -0.02630 -0.02587 -0.00969 -0.01747 0.47064

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/temp_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
0.04187 0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.04187 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.04187 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.02093 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.02093 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000 0.02093

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/born_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
3.35855 1.86892 1.87139 0.00233 0.00218 -0.00179
1.86892 3.37104 1.87285 0.00112 0.00085 -0.00007
1.87139 1.87285 3.37707 -0.00058 0.00038 -0.00057
0.00233 0.00112 -0.00058 1.88326 -0.00039 0.00065
0.00218 0.00085 0.00038 -0.00039 1.88229 0.00242
-0.00179 -0.00007 -0.00057 0.00065 0.00242 1.87968

118
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py Normal file
View File

@ -0,0 +1,118 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import sys
import numpy as np
def reduce_Born(Cf):
C = np.zeros((6,6), dtype=np.float64)
C[0,0] = Cf[0]
C[1,1] = Cf[1]
C[2,2] = Cf[2]
C[3,3] = Cf[3]
C[4,4] = Cf[4]
C[5,5] = Cf[5]
C[0,1] = Cf[6]
C[0,2] = Cf[7]
C[0,3] = Cf[8]
C[0,4] = Cf[9]
C[0,5] = Cf[10]
C[1,2] = Cf[11]
C[1,3] = Cf[12]
C[1,4] = Cf[13]
C[1,5] = Cf[14]
C[2,3] = Cf[15]
C[2,4] = Cf[16]
C[2,5] = Cf[17]
C[3,4] = Cf[18]
C[3,5] = Cf[19]
C[4,5] = Cf[20]
C = np.where(C,C,C.T)
return C
def compute_delta():
D = np.zeros((3,3,3,3))
for a in range(3):
for b in range(3):
for m in range(3):
for n in range(3):
D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
d = np.zeros((6,6))
d[0,0] = D[0,0,0,0]
d[1,1] = D[1,1,1,1]
d[2,2] = D[2,2,2,2]
d[3,3] = D[1,2,1,2]
d[4,4] = D[0,2,0,2]
d[5,5] = D[0,1,0,1]
d[0,1] = D[0,0,1,1]
d[0,2] = D[0,0,2,2]
d[0,3] = D[0,0,1,2]
d[0,4] = D[0,0,0,2]
d[0,5] = D[0,0,0,1]
d[1,2] = D[1,1,2,2]
d[1,3] = D[1,1,1,2]
d[1,4] = D[1,1,0,2]
d[1,5] = D[1,1,0,1]
d[2,3] = D[2,2,1,2]
d[2,4] = D[2,2,0,2]
d[2,5] = D[2,2,0,1]
d[3,4] = D[1,2,0,2]
d[3,5] = D[1,2,0,1]
d[4,5] = D[0,2,0,1]
d = np.where(d,d,d.T)
return d
def write_matrix(C, filename):
with open(filename, 'w') as f:
f.write("# Cij Matrix from post process computation\n")
for i in C:
f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
)
)
return
def main():
N = 500
vol = 27.047271**3 * 10**-30 # m^3
T = 60 # K
kb = 1.380649 * 10**-23 # J/K
kbT = T*kb # J
kcalmol2J = 4183.9954/(6.022*10**23)
born = np.loadtxt('born.out')
stre = np.loadtxt('vir.out')
stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa
try:
mean_born = np.mean(born[:, 1:], axis=0)
except IndexError:
mean_born = born[1:]
CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa
Cs = vol/kbT*np.cov(stre[:,1:].T)
Ct = N*kbT/vol * compute_delta()
C = CB - Cs + Ct
write_matrix(CB, 'born_matrix.out')
write_matrix(Cs, 'stre_matrix.out')
write_matrix(Ct, 'temp_matrix.out')
write_matrix(C, 'full_matrix.out')
C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
print(C*10**-9)
print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
return
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
raise SystemExit("User interruption.")

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/full_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
2.29922 1.17439 1.20854 -0.01653 -0.01684 0.01188
1.17439 2.20673 1.21718 -0.00781 -0.00753 0.00867
1.20854 1.21718 2.30804 0.01535 -0.01596 0.00426
-0.01653 -0.00781 0.01535 1.47647 -0.01355 -0.01601
-0.01684 -0.00753 -0.01596 -0.01355 1.37905 0.01975
0.01188 0.00867 0.00426 -0.01601 0.01975 1.40170

154
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov Normal file
View File

@ -0,0 +1,154 @@
# Numerical difference calculation
# of Born matrix
# Note that because of cubic symmetry and central forces, we have:
# C11, pure axial == positive mean value: 1,2,3
# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
# adjustable parameters
units real
variable nsteps index 100000 # length of run
variable nthermo index 1000 # thermo output interval
variable nlat equal 5 # size of box
variable T equal 60. # Temperature in K
variable rho equal 5.405 # Lattice spacing in A
atom_style atomic
lattice fcc ${rho}
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
create_box 1 box
create_atoms 1 box
mass * 39.948
velocity all create ${T} 87287 loop geom
velocity all zero linear
pair_style lj/cut 12.0
pair_coeff 1 1 0.238067 3.405
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
variable vol equal vol
thermo 100
fix aL all ave/time 1 1 1 v_vol ave running
fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
run 20000
unfix NPT
variable newL equal "f_aL^(1./3.)"
change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
unfix aL
reset_timestep 0
# Conversion variables
variable kb equal 1.38065e-23 # J/K
variable Myvol equal "vol*10^-30" # Volume in m^3
variable kbt equal "v_kb*v_T"
variable Nat equal atoms
variable Rhokbt equal "v_kbt*v_Nat/v_Myvol"
variable at2Pa equal 101325
variable kcalmol2J equal "4183.9954/(6.022e23)"
variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms
variable Pa2GPa equal 1e-9
# Born compute giving <C^b> terms
# The six virial stress component to compute <C^fl>
compute VIR all pressure NULL virial
compute born all born/matrix numdiff 1e-6 VIR
variable s1 equal "-c_VIR[1]*v_at2Pa"
variable s2 equal "-c_VIR[2]*v_at2Pa"
variable s3 equal "-c_VIR[3]*v_at2Pa"
variable s6 equal "-c_VIR[4]*v_at2Pa"
variable s5 equal "-c_VIR[5]*v_at2Pa"
variable s4 equal "-c_VIR[6]*v_at2Pa"
variable press equal press
# Average of Born term and vector to store stress
# for post processing
fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
fix APR all ave/time 1 1 1 v_press ave running
fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
thermo ${nthermo}
thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
thermo_modify line multi
fix 1 all nvt temp $T $T 100
run ${nsteps}
# Compute vector averages
# Note the indice switch.
# LAMMPS convention is NOT the Voigt notation.
variable aves1 equal "ave(f_VEC[1])"
variable aves2 equal "ave(f_VEC[2])"
variable aves3 equal "ave(f_VEC[3])"
variable aves4 equal "ave(f_VEC[6])"
variable aves5 equal "ave(f_VEC[5])"
variable aves6 equal "ave(f_VEC[4])"
# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
# is numerically instable. Here we go through the <(s-<s>)^2>
# definition.
# Computing difference relative to average values
variable ds1 vector "f_VEC[1]-v_aves1"
variable ds2 vector "f_VEC[2]-v_aves2"
variable ds3 vector "f_VEC[3]-v_aves3"
variable ds4 vector "f_VEC[4]-v_aves4"
variable ds5 vector "f_VEC[5]-v_aves5"
variable ds6 vector "f_VEC[6]-v_aves6"
# Squaring and averaging
variable dds1 vector "v_ds1*v_ds1"
variable dds2 vector "v_ds2*v_ds2"
variable dds3 vector "v_ds3*v_ds3"
variable vars1 equal "ave(v_dds1)"
variable vars2 equal "ave(v_dds2)"
variable vars3 equal "ave(v_dds3)"
variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
variable dds12 vector "v_ds1*v_ds2"
variable dds13 vector "v_ds1*v_ds3"
variable dds23 vector "v_ds2*v_ds3"
variable vars12 equal "ave(v_dds12)"
variable vars13 equal "ave(v_dds13)"
variable vars23 equal "ave(v_dds23)"
variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
variable dds4 vector "v_ds4*v_ds4"
variable dds5 vector "v_ds5*v_ds5"
variable dds6 vector "v_ds6*v_ds6"
variable vars4 equal "ave(v_dds4)"
variable vars5 equal "ave(v_dds5)"
variable vars6 equal "ave(v_dds6)"
variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
print """
C11 = ${aC11}
C12 = ${aC12}
C44 = ${aC44}
"""

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/stre_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
1.10120 0.69454 0.66285 0.01885 0.01902 -0.01367
0.69454 1.20617 0.65567 0.00893 0.00839 -0.00873
0.66285 0.65567 1.11090 -0.01593 0.01634 -0.00483
0.01885 0.00893 -0.01593 0.42772 0.01316 0.01666
0.01902 0.00839 0.01634 0.01316 0.52416 -0.01733
-0.01367 -0.00873 -0.00483 0.01666 -0.01733 0.49891

7
examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/temp_matrix.out Normal file
View File

@ -0,0 +1,7 @@
# Cij Matrix from post process computation
0.04187 0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 0.04187 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000 0.04187 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.02093 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000 0.02093 0.00000
0.00000 0.00000 0.00000 0.00000 0.00000 0.02093

13
examples/ELASTIC_T/BORN_MATRIX/Argon/README.md Normal file
View File

@ -0,0 +1,13 @@
This repository is a test case for the compute born/matrix. It provides short
scripts creating argon fcc crystal and computing the Born term using the two
methods described in the documentation.
In the __Analytical__ directory the terms are computed using the analytical
derivation of the Born term for the lj/cut pair style.
In the __Numdiff__ directory, the Born term is evaluated through small
numerical differences of the stress tensor. This method can be used with any
interaction potential.
Both script show examples on how to compute the full Cij elastic stiffness
tensor in LAMMPS.

1
examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/Si.sw

68
examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in Normal file
View File

@ -0,0 +1,68 @@
# Average moduli for cubic crystals
variable C11cubic equal (${C11}+${C22}+${C33})/3.0
variable C12cubic equal (${C12}+${C13}+${C23})/3.0
variable C44cubic equal (${C44}+${C55}+${C66})/3.0
variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0
variable shearmodulus1 equal ${C44cubic}
variable shearmodulus2 equal (${C11cubic}-${C12cubic})/2.0
variable poissonratio equal 1.0/(1.0+${C11cubic}/${C12cubic})
# For Stillinger-Weber silicon, the analytical results
# are known to be (E. R. Cowley, 1988):
# C11 = 151.4 GPa
# C12 = 76.4 GPa
# C44 = 56.4 GPa
#print "========================================="
#print "Components of the Elastic Constant Tensor"
#print "========================================="
print "Elastic Constant C11 = ${C11} ${cunits}"
print "Elastic Constant C22 = ${C22} ${cunits}"
print "Elastic Constant C33 = ${C33} ${cunits}"
print "Elastic Constant C12 = ${C12} ${cunits}"
print "Elastic Constant C13 = ${C13} ${cunits}"
print "Elastic Constant C23 = ${C23} ${cunits}"
print "Elastic Constant C44 = ${C44} ${cunits}"
print "Elastic Constant C55 = ${C55} ${cunits}"
print "Elastic Constant C66 = ${C66} ${cunits}"
print "Elastic Constant C14 = ${C14} ${cunits}"
print "Elastic Constant C15 = ${C15} ${cunits}"
print "Elastic Constant C16 = ${C16} ${cunits}"
print "Elastic Constant C24 = ${C24} ${cunits}"
print "Elastic Constant C25 = ${C25} ${cunits}"
print "Elastic Constant C26 = ${C26} ${cunits}"
print "Elastic Constant C34 = ${C34} ${cunits}"
print "Elastic Constant C35 = ${C35} ${cunits}"
print "Elastic Constant C36 = ${C36} ${cunits}"
print "Elastic Constant C45 = ${C45} ${cunits}"
print "Elastic Constant C46 = ${C46} ${cunits}"
print "Elastic Constant C56 = ${C56} ${cunits}"
print "========================================="
print "Average properties for a cubic crystal"
print "========================================="
print "Bulk Modulus = ${bulkmodulus} ${cunits}"
print "Shear Modulus 1 = ${shearmodulus1} ${cunits}"
print "Shear Modulus 2 = ${shearmodulus2} ${cunits}"
print "Poisson Ratio = ${poissonratio}"
# summarize sampling protocol
variable tmp equal atoms
print "Number of atoms = ${tmp}"
print "Stress sampling interval = ${nevery}"
variable tmp equal ${nrun}/${nevery}
print "Stress sample count = ${tmp}"
print "Born sampling interval = ${neveryborn}"
variable tmp equal ${nrun}/${neveryborn}
print "Born sample count = ${tmp}"

25
examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic Normal file
View File

@ -0,0 +1,25 @@
# Compute elastic constant tensor for a crystal at finite temperature
# These settings replicate the 1477~K benchmark of
# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
# here: Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
include init.in
# Compute initial state
include potential.in
thermo_style custom step temp pe press density
run ${nequil}
# Run dynamics
include potential.in
include output.in
run ${nrun}
# Output final values
include final_output.in

59
examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in Normal file
View File

@ -0,0 +1,59 @@
# NOTE: This script can be modified for different atomic structures,
# units, etc. See in.elastic for more info.
#
# Define MD parameters
# These can be modified by the user
# These settings replicate the 1477~K benchmark of
# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
# select temperature and pressure (lattice constant)
variable temp index 1477.0 # temperature of initial sample
variable a index 5.457 # lattice constant
# select sampling parameters, important for speed/convergence
variable nthermo index 1500 # interval for thermo output
variable nevery index 10 # stress sampling interval
variable neveryborn index 100 # Born sampling interval
variable timestep index 0.000766 # timestep
variable nlat index 3 # number of lattice unit cells
# other settings
variable mass1 index 28.06 # mass
variable tdamp index 0.01 # time constant for thermostat
variable seed index 123457 # seed for thermostat
variable thermostat index 1 # 0 if NVE, 1 if NVT
variable delta index 1.0e-6 # Born numdiff strain magnitude
# hard-coded rules-of-thumb for run length, etc.
variable nfreq equal ${nthermo} # interval for averaging output
variable nrepeat equal floor(${nfreq}/${nevery}) # number of samples
variable nrepeatborn equal floor(${nfreq}/${neveryborn}) # number of samples
variable nequil equal 10*${nthermo} # length of equilibration run
variable nrun equal 100*${nthermo} # length of equilibrated run
# generate the box and atom positions using a diamond lattice
units metal
boundary p p p
# this generates a standard 8-atom cubic cell
lattice diamond $a
region box prism 0 1 0 1 0 1 0 0 0
# this generates a 2-atom triclinic cell
#include tri.in
create_box 1 box
create_atoms 1 box
mass 1 ${mass1}
replicate ${nlat} ${nlat} ${nlat}
velocity all create ${temp} 87287

140
examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in Normal file
View File

@ -0,0 +1,140 @@
# Setup output
# Stress fluctuation term F
compute stress all pressure thermo_temp
variable s1 equal c_stress[1]
variable s2 equal c_stress[2]
variable s3 equal c_stress[3]
variable s4 equal c_stress[6]
variable s5 equal c_stress[5]
variable s6 equal c_stress[4]
variable s11 equal v_s1*v_s1
variable s22 equal v_s2*v_s2
variable s33 equal v_s3*v_s3
variable s44 equal v_s4*v_s4
variable s55 equal v_s5*v_s5
variable s66 equal v_s6*v_s6
variable s33 equal v_s3*v_s3
variable s12 equal v_s1*v_s2
variable s13 equal v_s1*v_s3
variable s14 equal v_s1*v_s4
variable s15 equal v_s1*v_s5
variable s16 equal v_s1*v_s6
variable s23 equal v_s2*v_s3
variable s24 equal v_s2*v_s4
variable s25 equal v_s2*v_s5
variable s26 equal v_s2*v_s6
variable s34 equal v_s3*v_s4
variable s35 equal v_s3*v_s5
variable s36 equal v_s3*v_s6
variable s45 equal v_s4*v_s5
variable s46 equal v_s4*v_s6
variable s56 equal v_s5*v_s6
variable mytemp equal temp
variable mypress equal press
variable mype equal pe/atoms
fix avt all ave/time ${nevery} ${nrepeat} ${nfreq} v_mytemp ave running
fix avp all ave/time ${nevery} ${nrepeat} ${nfreq} v_mypress ave running
fix avpe all ave/time ${nevery} ${nrepeat} ${nfreq} v_mype ave running
fix avs all ave/time ${nevery} ${nrepeat} ${nfreq} v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 ave running
fix avssq all ave/time ${nevery} ${nrepeat} ${nfreq} &
v_s11 v_s22 v_s33 v_s44 v_s55 v_s66 &
v_s12 v_s13 v_s14 v_s15 v_s16 &
v_s23 v_s24 v_s25 v_s26 &
v_s34 v_s35 v_s36 &
v_s45 v_s46 &
v_s56 &
ave running
# bar to GPa
variable pconv equal 1.0e5/1.0e9
variable cunits index GPa
# metal unit constants from LAMMPS
# force->nktv2p = 1.6021765e6;
# force->boltz = 8.617343e-5;
variable boltz equal 8.617343e-5
variable nktv2p equal 1.6021765e6
variable vkt equal vol/(${boltz}*${temp})/${nktv2p}
variable ffac equal ${pconv}*${vkt}
variable F11 equal -(f_avssq[1]-f_avs[1]*f_avs[1])*${ffac}
variable F22 equal -(f_avssq[2]-f_avs[2]*f_avs[2])*${ffac}
variable F33 equal -(f_avssq[3]-f_avs[3]*f_avs[3])*${ffac}
variable F44 equal -(f_avssq[4]-f_avs[4]*f_avs[4])*${ffac}
variable F55 equal -(f_avssq[5]-f_avs[5]*f_avs[5])*${ffac}
variable F66 equal -(f_avssq[6]-f_avs[6]*f_avs[6])*${ffac}
variable F12 equal -(f_avssq[7]-f_avs[1]*f_avs[2])*${ffac}
variable F13 equal -(f_avssq[8]-f_avs[1]*f_avs[3])*${ffac}
variable F14 equal -(f_avssq[9]-f_avs[1]*f_avs[4])*${ffac}
variable F15 equal -(f_avssq[10]-f_avs[1]*f_avs[5])*${ffac}
variable F16 equal -(f_avssq[11]-f_avs[1]*f_avs[6])*${ffac}
variable F23 equal -(f_avssq[12]-f_avs[2]*f_avs[3])*${ffac}
variable F24 equal -(f_avssq[13]-f_avs[2]*f_avs[4])*${ffac}
variable F25 equal -(f_avssq[14]-f_avs[2]*f_avs[5])*${ffac}
variable F26 equal -(f_avssq[15]-f_avs[2]*f_avs[6])*${ffac}
variable F34 equal -(f_avssq[16]-f_avs[3]*f_avs[4])*${ffac}
variable F35 equal -(f_avssq[17]-f_avs[3]*f_avs[5])*${ffac}
variable F36 equal -(f_avssq[18]-f_avs[3]*f_avs[6])*${ffac}
variable F45 equal -(f_avssq[19]-f_avs[4]*f_avs[5])*${ffac}
variable F46 equal -(f_avssq[20]-f_avs[4]*f_avs[6])*${ffac}
variable F56 equal -(f_avssq[21]-f_avs[5]*f_avs[6])*${ffac}
# Born term
compute virial all pressure NULL virial
compute born all born/matrix numdiff ${delta} virial
fix avborn all ave/time ${neveryborn} ${nrepeatborn} ${nfreq} c_born[*] ave running
variable bfac equal ${pconv}*${nktv2p}/vol
variable B vector f_avborn*${bfac}
# Kinetic term
variable kfac equal ${pconv}*${nktv2p}*atoms*${boltz}*${temp}/vol
variable K11 equal 4.0*${kfac}
variable K22 equal 4.0*${kfac}
variable K33 equal 4.0*${kfac}
variable K44 equal 2.0*${kfac}
variable K55 equal 2.0*${kfac}
variable K66 equal 2.0*${kfac}
# Add F, K, and B together
variable C11 equal v_F11+v_B[1]+v_K11
variable C22 equal v_F22+v_B[2]+v_K22
variable C33 equal v_F33+v_B[3]+v_K33
variable C44 equal v_F44+v_B[4]+v_K44
variable C55 equal v_F55+v_B[5]+v_K55
variable C66 equal v_F66+v_B[6]+v_K66
variable C12 equal v_F12+v_B[7]
variable C13 equal v_F13+v_B[8]
variable C14 equal v_F14+v_B[9]
variable C15 equal v_F15+v_B[10]
variable C16 equal v_F16+v_B[11]
variable C23 equal v_F23+v_B[12]
variable C24 equal v_F24+v_B[13]
variable C25 equal v_F25+v_B[14]
variable C26 equal v_F26+v_B[15]
variable C34 equal v_F34+v_B[16]
variable C35 equal v_F35+v_B[17]
variable C36 equal v_F36+v_B[18]
variable C45 equal v_F45+v_B[19]
variable C46 equal v_F46+v_B[20]
variable C56 equal v_F56+v_B[21]
thermo ${nthermo}
thermo_style custom step temp pe press density f_avt f_avp f_avpe v_F11 v_F22 v_F33 v_F44 v_F55 v_F66 v_F12 v_F13 v_F23 v_B[1] v_B[2] v_B[3] v_B[4] v_B[5] v_B[6] v_B[7] v_B[8] v_B[12]
thermo_modify norm no

21
examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in Normal file
View File

@ -0,0 +1,21 @@
# NOTE: This script can be modified for different pair styles
# See in.elastic for more info.
reset_timestep 0
# Choose potential
pair_style sw
pair_coeff * * Si.sw Si
# Setup neighbor style
neighbor 1.0 nsq
neigh_modify once no every 1 delay 0 check yes
# Setup MD
timestep ${timestep}
fix 4 all nve
if "${thermostat} == 1" then &
"fix 5 all langevin ${temp} ${temp} ${tdamp} ${seed}"

26
examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in Normal file
View File

@ -0,0 +1,26 @@
# this generates a 2-atom triclinic cell
# due to rotation on to x-axis,
# elastic constant analysis is not working yet
# unit lattice vectors are
# a1 = (1 0 0)
# a2 = (1/2 sqrt3/2 0)
# a3 = (1/2 1/(2sqrt3) sqrt2/sqrt3)
variable a1x equal 1
variable a2x equal 1/2
variable a2y equal sqrt(3)/2
variable a3x equal 1/2
variable a3y equal 1/(2*sqrt(3))
variable a3z equal sqrt(2/3)
variable l equal $a/sqrt(2)
lattice custom ${l} &
a1 ${a1x} 0 0 &
a2 ${a2x} ${a2y} 0.0 &
a3 ${a3x} ${a3y} ${a3z} &
basis 0 0 0 &
basis 0.25 0.25 0.25 &
spacing 1 1 1
region box prism 0 ${a1x} 0 ${a2y} 0 ${a3z} ${a2x} ${a3x} ${a3y}

1
examples/ELASTIC_T/DEFORMATION/Silicon/Si.sw Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/Si.sw

0
examples/ELASTIC_T/displace.mod → examples/ELASTIC_T/DEFORMATION/Silicon/displace.mod
View File

0
examples/ELASTIC_T/in.elastic → examples/ELASTIC_T/DEFORMATION/Silicon/in.elastic
View File

0
examples/ELASTIC_T/init.mod → examples/ELASTIC_T/DEFORMATION/Silicon/init.mod
View File

2
examples/ELASTIC_T/potential.mod → examples/ELASTIC_T/DEFORMATION/Silicon/potential.mod
View File

@ -2,7 +2,7 @@
# See in.elastic for more info.
# we must undefine any fix ave/* fix before using reset_timestep
if "$(is_defined(fix,avp)" then "unfix avp"
if "$(is_defined(fix,avp))" then "unfix avp"
reset_timestep 0
# Choose potential

18
examples/ELASTIC_T/Si.sw
View File

@ -1,18 +0,0 @@
# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber, Phys Rev B, 31, 5262, (1985)
# Stillinger-Weber parameters for various elements and mixtures
# multiple entries can be added to this file, LAMMPS reads the ones it needs
# these entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
# format of a single entry (one or more lines):
# element 1, element 2, element 3,
# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
# Here are the original parameters in metal units, for Silicon from:
#
# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985)
#
Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333
7.049556277 0.6022245584 4.0 0.0 0.0

17
examples/README
View File

@ -94,7 +94,7 @@ msst: MSST shock dynamics
nb3b: use of nonbonded 3-body harmonic pair style
neb: nudged elastic band (NEB) calculation for barrier finding
nemd: non-equilibrium MD of 2d sheared system
numdiff: numerical difference computation of forces and virial
numdiff: numerical difference computation of forces, virial, and Born matrix
obstacle: flow around two voids in a 2d channel
peptide: dynamics of a small solvated peptide chain (5-mer)
peri: Peridynamic model of cylinder impacted by indenter
@ -153,12 +153,19 @@ either by itself or in tandem with another code or library. See the
COUPLE/README file to get started.
The ELASTIC directory has an example script for computing elastic
constants at zero temperature, using an Si example. See the
stiffness tensor (elastic constants)
at zero temperature, using an Si example. See the
ELASTIC/in.elastic file for more info.
The ELASTIC_T directory has an example script for computing elastic
constants at finite temperature, using an Si example. See the
ELASTIC_T/in.elastic file for more info.
The ELASTIC_T directory has example scripts for the computing elastic
stiffness tensor at finite temperature. Two different methods are
demonstrated. DEFORMATION estimates the change in the average
stress tensor between multiple simulations
in which small finite deformations are made to the simulation cell.
BORN_MATRIX runs a single simulation in which the Born matrix and stress
fluctuations are averaged. The second method
is newer in LAMMPS and is generally more efficient and
more reliable.
The HEAT directory has example scripts for heat exchange algorithms
(e.g. used for establishing a thermal gradient), using two different

11
examples/numdiff/README.md
View File

@ -1,8 +1,11 @@
# LAMMPS FIX NUMDIFF EXAMPLE
# LAMMPS NUMDIFF EXAMPLES FOR FORCES, VIRIAL, and BORN MATRIX
## Numerical Difference Fix
## Numerical Difference Fixes and Computes
This directory contains the ingredients to run an NVE simulation using the numerical difference fix and calculate error in forces.
This directory contains the input script for an NVE simulation with
fix numdiff, fix numdiff/virial, and compute born/matrix numdiff.
In each cases, results are compared to exact analytic expressions
and the small relative differences are reported.
Example:
```
@ -10,4 +13,4 @@ NP=4 #number of processors
mpirun -np $NP lmp_mpi -in.numdiff
```
## Required LAMMPS packages: MOLECULE package
## Required LAMMPS packages: EXTRA-FIX, EXTRA-COMPUTE

67
examples/numdiff/in.numdiff
View File

@ -1,5 +1,5 @@
# Numerical difference calculation
# of error in forces and virial stress
# of error in forces, virial stress, and Born matrix
# adjustable parameters
@ -8,8 +8,10 @@ variable nthermo index 10 # thermo output interval
variable ndump index 500 # dump output interval
variable nlat index 3 # size of box
variable fdelta index 1.0e-4 # displacement size
variable vdelta index 1.0e-6 # strain size
variable vdelta index 1.0e-6 # strain size for numdiff/virial
variable bdelta index 1.0e-8 # strain size for numdiff Born matrix
variable temp index 10.0 # temperature
variable nugget equal 1.0e-6 # regularization for relerr
units metal
atom_style atomic
@ -23,8 +25,8 @@ mass 1 39.903
velocity all create ${temp} 2357 mom yes dist gaussian
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
pair_style lj/cut 5.0
pair_coeff 1 1 0.0102701 3.42
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
@ -42,7 +44,7 @@ variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
compute faverrsq all reduce ave v_ferrsq
variable fsq atom fx^2+fy^2+fz^2
compute favsq all reduce ave v_fsq
variable frelerr equal sqrt(c_faverrsq/c_favsq)
variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget}))
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
# define numerical virial stress tensor calculation
@ -67,8 +69,59 @@ variable vsq equal "c_myvirial[1]^2 + &
c_myvirial[4]^2 + &
c_myvirial[5]^2 + &
c_myvirial[6]^2"
variable vrelerr equal sqrt(v_verrsq/v_vsq)
variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget}))
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr
# define numerical Born matrix calculation
compute bornnum all born/matrix numdiff ${bdelta} myvirial
compute born all born/matrix
variable berr vector c_bornnum-c_born
variable berrsq equal "v_berr[1]^2 + &
v_berr[2]^2 + &
v_berr[3]^2 + &
v_berr[4]^2 + &
v_berr[5]^2 + &
v_berr[6]^2 + &
v_berr[7]^2 + &
v_berr[8]^2 + &
v_berr[9]^2 + &
v_berr[10]^2 + &
v_berr[11]^2 + &
v_berr[12]^2 + &
v_berr[13]^2 + &
v_berr[14]^2 + &
v_berr[15]^2 + &
v_berr[16]^2 + &
v_berr[17]^2 + &
v_berr[18]^2 + &
v_berr[19]^2 + &
v_berr[20]^2 + &
v_berr[21]^2"
variable bsq equal "c_born[1]^2 + &
c_born[2]^2 + &
c_born[3]^2 + &
c_born[4]^2 + &
c_born[5]^2 + &
c_born[6]^2 + &
c_born[7]^2 + &
c_born[8]^2 + &
c_born[9]^2 + &
c_born[10]^2 + &
c_born[11]^2 + &
c_born[12]^2 + &
c_born[13]^2 + &
c_born[14]^2 + &
c_born[15]^2 + &
c_born[16]^2 + &
c_born[17]^2 + &
c_born[18]^2 + &
c_born[19]^2 + &
c_born[20]^2 + &
c_born[21]^2"
variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget}))
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr
thermo ${nthermo}
run ${nsteps}

197
examples/numdiff/log.19Feb2022.log.numdiff.g++.1 Normal file
View File

@ -0,0 +1,197 @@
LAMMPS (17 Feb 2022)
# Numerical difference calculation
# of error in forces, virial stress, and Born matrix
# adjustable parameters
variable nsteps index 500 # length of run
variable nthermo index 10 # thermo output interval
variable ndump index 500 # dump output interval
variable nlat index 3 # size of box
variable fdelta index 1.0e-4 # displacement size
variable vdelta index 1.0e-6 # strain size for numdiff/virial
variable bdelta index 1.0e-8 # strain size for numdiff Born matrix
variable temp index 10.0 # temperature
variable nugget equal 1.0e-6 # regularization for relerr
units metal
atom_style atomic
atom_modify map yes
lattice fcc 5.358000
Lattice spacing in x,y,z = 5.358 5.358 5.358
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
region box block 0 3 0 ${nlat} 0 ${nlat}
region box block 0 3 0 3 0 ${nlat}
region box block 0 3 0 3 0 3
create_box 1 box
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 108 atoms
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
create_atoms CPU = 0.000 seconds
mass 1 39.903
velocity all create ${temp} 2357 mom yes dist gaussian
velocity all create 10.0 2357 mom yes dist gaussian
pair_style lj/cut 5.0
pair_coeff 1 1 0.0102701 3.42
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
timestep 0.001
fix nve all nve
# define numerical force calculation
fix numforce all numdiff ${nthermo} ${fdelta}
fix numforce all numdiff 10 ${fdelta}
fix numforce all numdiff 10 1.0e-4
variable ferrx atom f_numforce[1]-fx
variable ferry atom f_numforce[2]-fy
variable ferrz atom f_numforce[3]-fz
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
compute faverrsq all reduce ave v_ferrsq
variable fsq atom fx^2+fy^2+fz^2
compute favsq all reduce ave v_fsq
variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget}))
variable frelerr equal sqrt(c_faverrsq/(c_favsq+1e-06))
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
# define numerical virial stress tensor calculation
compute myvirial all pressure NULL virial
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
fix numvirial all numdiff/virial 10 ${vdelta}
fix numvirial all numdiff/virial 10 1.0e-6
variable errxx equal f_numvirial[1]-c_myvirial[1]
variable erryy equal f_numvirial[2]-c_myvirial[2]
variable errzz equal f_numvirial[3]-c_myvirial[3]
variable erryz equal f_numvirial[4]-c_myvirial[6]
variable errxz equal f_numvirial[5]-c_myvirial[5]
variable errxy equal f_numvirial[6]-c_myvirial[4]
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget}))
variable vrelerr equal sqrt(v_verrsq/(v_vsq+1e-06))
# define numerical Born matrix calculation
compute bornnum all born/matrix numdiff ${bdelta} myvirial
compute bornnum all born/matrix numdiff 1.0e-8 myvirial
compute born all born/matrix
variable berr vector c_bornnum-c_born
variable berrsq equal "v_berr[1]^2 + v_berr[2]^2 + v_berr[3]^2 + v_berr[4]^2 + v_berr[5]^2 + v_berr[6]^2 + v_berr[7]^2 + v_berr[8]^2 + v_berr[9]^2 + v_berr[10]^2 + v_berr[11]^2 + v_berr[12]^2 + v_berr[13]^2 + v_berr[14]^2 + v_berr[15]^2 + v_berr[16]^2 + v_berr[17]^2 + v_berr[18]^2 + v_berr[19]^2 + v_berr[20]^2 + v_berr[21]^2"
variable bsq equal "c_born[1]^2 + c_born[2]^2 + c_born[3]^2 + c_born[4]^2 + c_born[5]^2 + c_born[6]^2 + c_born[7]^2 + c_born[8]^2 + c_born[9]^2 + c_born[10]^2 + c_born[11]^2 + c_born[12]^2 + c_born[13]^2 + c_born[14]^2 + c_born[15]^2 + c_born[16]^2 + c_born[17]^2 + c_born[18]^2 + c_born[19]^2 + c_born[20]^2 + c_born[21]^2"
variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget}))
variable brelerr equal sqrt(v_berrsq/(v_bsq+1e-06))
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr
thermo ${nthermo}
thermo 10
run ${nsteps}
run 500
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5
ghost atom cutoff = 5
binsize = 2.5, bins = 7 7 7
2 neighbor lists, perpetual/occasional/extra = 1 1 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
(2) compute born/matrix, occasional, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 6.823 | 6.823 | 6.823 Mbytes
Step Temp PotEng TotEng Press v_frelerr v_vrelerr v_brelerr
0 10 -6.6101864 -6.471878 947.70558 5.7012262e-09 1.4849734e-08 9.036398e-09
10 9.9357441 -6.6092976 -6.471878 949.3486 1.3828998e-08 1.9248385e-09 4.0233493e-09
20 9.7444561 -6.6066519 -6.4718779 954.23637 1.385204e-08 1.7476399e-09 4.0061081e-09
30 9.4311148 -6.6023181 -6.4718779 962.23331 1.4147226e-08 1.7647816e-09 3.3866543e-09
40 9.0043293 -6.5964152 -6.4718778 973.10762 1.4128155e-08 1.6390138e-09 3.2652821e-09
50 8.4762135 -6.5891108 -6.4718777 986.53572 1.4168048e-08 2.3910821e-09 4.7266627e-09
60 7.8621735 -6.5806179 -6.4718775 1002.1092 1.411958e-08 2.0683414e-09 2.6951001e-09
70 7.1805874 -6.5711908 -6.4718773 1019.3448 1.4139911e-08 1.6084571e-09 3.1477301e-09
80 6.4523557 -6.5611186 -6.4718771 1037.6974 1.4105096e-08 1.9929271e-09 3.4733802e-09
90 5.7003071 -6.5507169 -6.4718769 1056.5767 1.4084183e-08 1.750579e-09 4.310104e-09
100 4.9484503 -6.5403179 -6.4718767 1075.3674 1.4063796e-08 1.0250271e-09 2.9213594e-09
110 4.221081 -6.5302576 -6.4718765 1093.4526 1.400901e-08 1.389277e-09 4.3909721e-09
120 3.5417733 -6.520862 -6.4718763 1110.2388 1.4038158e-08 8.6231891e-10 2.5890696e-09
130 2.9323072 -6.5124324 -6.4718762 1125.183 1.4048645e-08 7.0840985e-10 3.388192e-09
140 2.411607 -6.5052306 -6.471876 1137.8182 1.3968429e-08 1.8508015e-09 3.2976031e-09
150 1.9947801 -6.4994654 -6.4718759 1147.7764 1.395965e-08 1.9484728e-09 4.2924605e-09
160 1.6923481 -6.4952825 -6.4718759 1154.8063 1.3948606e-08 1.5275137e-09 4.0204309e-09
170 1.5097515 -6.492757 -6.4718759 1158.7853 1.3845523e-08 1.5455e-09 4.8781309e-09
180 1.4471795 -6.4918916 -6.4718759 1159.7221 1.3788451e-08 1.578099e-09 3.0795316e-09
190 1.4997431 -6.4926187 -6.471876 1157.7529 1.374841e-08 2.142073e-09 2.4376961e-09
200 1.6579637 -6.4948072 -6.4718761 1153.1286 1.3674788e-08 2.111894e-09 3.7055708e-09
210 1.908522 -6.4982727 -6.4718763 1146.1965 1.3639408e-08 1.2386489e-09 3.160881e-09
220 2.23518 -6.5027908 -6.4718764 1137.3775 1.3524209e-08 1.7016573e-09 3.6982265e-09
230 2.6197892 -6.5081105 -6.4718766 1127.1415 1.3344007e-08 1.5843477e-09 3.7272821e-09
240 3.043298 -6.5139681 -6.4718768 1115.9815 1.3245227e-08 1.5502368e-09 3.898015e-09
250 3.4866901 -6.5201007 -6.4718769 1104.3906 1.3080142e-08 1.369987e-09 4.9133863e-09
260 3.9318061 -6.5262572 -6.4718771 1092.84 1.2885339e-08 1.0743728e-09 5.7271364e-09
270 4.3620216 -6.5322076 -6.4718772 1081.7617 1.2705966e-08 1.3618619e-09 2.3225062e-09
280 4.7627723 -6.5377504 -6.4718773 1071.5341 1.2480463e-08 1.4346869e-09 3.281167e-09
290 5.1219322 -6.542718 -6.4718774 1062.4716 1.2434727e-08 2.1935942e-09 2.8198924e-09
300 5.4300557 -6.5469796 -6.4718774 1054.8177 1.2321314e-08 8.2365886e-10 3.2731015e-09
310 5.6804997 -6.5504435 -6.4718774 1048.7409 1.2300884e-08 1.4855741e-09 4.1031988e-09
320 5.8694423 -6.5530567 -6.4718774 1044.3341 1.2483087e-08 1.8711589e-09 3.9368436e-09
330 5.9958115 -6.5548045 -6.4718774 1041.6165 1.2627617e-08 1.9256986e-09 4.3283764e-09
340 6.0611353 -6.555708 -6.4718774 1040.5369 1.2935701e-08 1.6609255e-09 3.8728039e-09
350 6.0693222 -6.5558211 -6.4718773 1040.9803 1.3218179e-08 1.985355e-09 2.618577e-09
360 6.0263776 -6.5552271 -6.4718773 1042.7755 1.3471701e-08 1.5125203e-09 2.936238e-09
370 5.9400629 -6.5540332 -6.4718772 1045.7049 1.3676495e-08 1.7364093e-09 2.9097362e-09
380 5.8195019 -6.5523657 -6.4718771 1049.515 1.3859995e-08 1.6834835e-09 2.7416302e-09
390 5.6747442 -6.5503635 -6.471877 1053.9288 1.3987553e-08 1.7893896e-09 2.8552537e-09
400 5.5162948 -6.5481719 -6.4718769 1058.6583 1.4091878e-08 1.4468098e-09 3.2733654e-09
410 5.3546269 -6.5459358 -6.4718768 1063.4182 1.4188438e-08 1.7231047e-09 3.3165187e-09
420 5.1996958 -6.5437929 -6.4718768 1067.9384 1.4205207e-08 1.3551982e-09 3.8687611e-09
430 5.0604771 -6.5418673 -6.4718767 1071.9767 1.4267199e-08 1.361845e-09 3.1210672e-09
440 4.9445529 -6.5402639 -6.4718766 1075.3292 1.4253464e-08 1.3945282e-09 2.6483572e-09
450 4.8577717 -6.5390637 -6.4718766 1077.8394 1.4240998e-08 1.8767323e-09 3.2040422e-09
460 4.8040023 -6.53832 -6.4718766 1079.4048 1.4242259e-08 1.4785379e-09 3.4402279e-09
470 4.7849977 -6.5380571 -6.4718766 1079.9795 1.4227939e-08 1.8623848e-09 4.3634918e-09
480 4.8003794 -6.5382699 -6.4718766 1079.5756 1.4215836e-08 1.2821795e-09 2.6846581e-09
490 4.8477405 -6.538925 -6.4718767 1078.2596 1.4186541e-08 2.47604e-09 3.2044632e-09
500 4.9228588 -6.539964 -6.4718767 1076.1469 1.4099819e-08 1.6653302e-09 3.267113e-09
Loop time of 0.458483 on 1 procs for 500 steps with 108 atoms
Performance: 94.224 ns/day, 0.255 hours/ns, 1090.552 timesteps/s
99.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.0042278 | 0.0042278 | 0.0042278 | 0.0 | 0.92
Neigh | 0.02481 | 0.02481 | 0.02481 | 0.0 | 5.41
Comm | 0.002944 | 0.002944 | 0.002944 | 0.0 | 0.64
Output | 0.014731 | 0.014731 | 0.014731 | 0.0 | 3.21
Modify | 0.41122 | 0.41122 | 0.41122 | 0.0 | 89.69
Other | | 0.0005545 | | | 0.12
Nlocal: 108 ave 108 max 108 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 256 ave 256 max 256 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 648 ave 648 max 648 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 648
Ave neighs/atom = 6
Neighbor list builds = 500
Dangerous builds not checked
Total wall time: 0:00:00

197
examples/numdiff/log.19Feb2022.log.numdiff.g++.4 Normal file
View File

@ -0,0 +1,197 @@
LAMMPS (17 Feb 2022)
# Numerical difference calculation
# of error in forces, virial stress, and Born matrix
# adjustable parameters
variable nsteps index 500 # length of run
variable nthermo index 10 # thermo output interval
variable ndump index 500 # dump output interval
variable nlat index 3 # size of box
variable fdelta index 1.0e-4 # displacement size
variable vdelta index 1.0e-6 # strain size for numdiff/virial
variable bdelta index 1.0e-8 # strain size for numdiff Born matrix
variable temp index 10.0 # temperature
variable nugget equal 1.0e-6 # regularization for relerr
units metal
atom_style atomic
atom_modify map yes
lattice fcc 5.358000
Lattice spacing in x,y,z = 5.358 5.358 5.358
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
region box block 0 3 0 ${nlat} 0 ${nlat}
region box block 0 3 0 3 0 ${nlat}
region box block 0 3 0 3 0 3
create_box 1 box
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 108 atoms
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
create_atoms CPU = 0.000 seconds
mass 1 39.903
velocity all create ${temp} 2357 mom yes dist gaussian
velocity all create 10.0 2357 mom yes dist gaussian
pair_style lj/cut 5.0
pair_coeff 1 1 0.0102701 3.42
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
timestep 0.001
fix nve all nve
# define numerical force calculation
fix numforce all numdiff ${nthermo} ${fdelta}
fix numforce all numdiff 10 ${fdelta}
fix numforce all numdiff 10 1.0e-4
variable ferrx atom f_numforce[1]-fx
variable ferry atom f_numforce[2]-fy
variable ferrz atom f_numforce[3]-fz
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
compute faverrsq all reduce ave v_ferrsq
variable fsq atom fx^2+fy^2+fz^2
compute favsq all reduce ave v_fsq
variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget}))
variable frelerr equal sqrt(c_faverrsq/(c_favsq+1e-06))
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
# define numerical virial stress tensor calculation
compute myvirial all pressure NULL virial
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
fix numvirial all numdiff/virial 10 ${vdelta}
fix numvirial all numdiff/virial 10 1.0e-6
variable errxx equal f_numvirial[1]-c_myvirial[1]
variable erryy equal f_numvirial[2]-c_myvirial[2]
variable errzz equal f_numvirial[3]-c_myvirial[3]
variable erryz equal f_numvirial[4]-c_myvirial[6]
variable errxz equal f_numvirial[5]-c_myvirial[5]
variable errxy equal f_numvirial[6]-c_myvirial[4]
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget}))
variable vrelerr equal sqrt(v_verrsq/(v_vsq+1e-06))
# define numerical Born matrix calculation
compute bornnum all born/matrix numdiff ${bdelta} myvirial
compute bornnum all born/matrix numdiff 1.0e-8 myvirial
compute born all born/matrix
variable berr vector c_bornnum-c_born
variable berrsq equal "v_berr[1]^2 + v_berr[2]^2 + v_berr[3]^2 + v_berr[4]^2 + v_berr[5]^2 + v_berr[6]^2 + v_berr[7]^2 + v_berr[8]^2 + v_berr[9]^2 + v_berr[10]^2 + v_berr[11]^2 + v_berr[12]^2 + v_berr[13]^2 + v_berr[14]^2 + v_berr[15]^2 + v_berr[16]^2 + v_berr[17]^2 + v_berr[18]^2 + v_berr[19]^2 + v_berr[20]^2 + v_berr[21]^2"
variable bsq equal "c_born[1]^2 + c_born[2]^2 + c_born[3]^2 + c_born[4]^2 + c_born[5]^2 + c_born[6]^2 + c_born[7]^2 + c_born[8]^2 + c_born[9]^2 + c_born[10]^2 + c_born[11]^2 + c_born[12]^2 + c_born[13]^2 + c_born[14]^2 + c_born[15]^2 + c_born[16]^2 + c_born[17]^2 + c_born[18]^2 + c_born[19]^2 + c_born[20]^2 + c_born[21]^2"
variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget}))
variable brelerr equal sqrt(v_berrsq/(v_bsq+1e-06))
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr
thermo ${nthermo}
thermo 10
run ${nsteps}
run 500
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5
ghost atom cutoff = 5
binsize = 2.5, bins = 7 7 7
2 neighbor lists, perpetual/occasional/extra = 1 1 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
(2) compute born/matrix, occasional, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 6.816 | 6.816 | 6.816 Mbytes
Step Temp PotEng TotEng Press v_frelerr v_vrelerr v_brelerr
0 10 -6.6101864 -6.471878 947.70558 1.9110624e-09 9.4407596e-10 3.1867416e-09
10 9.9369961 -6.6093149 -6.471878 949.31222 1.3055176e-08 4.996456e-10 2.7421655e-09
20 9.7500224 -6.6067289 -6.4718779 954.07898 1.3721178e-08 5.6039795e-10 2.3689718e-09
30 9.4448115 -6.6025075 -6.4718779 961.85502 1.3813156e-08 6.8451692e-10 1.9844663e-09
40 9.0305392 -6.5967776 -6.4718777 972.39819 1.3961749e-08 3.1134064e-10 1.7915052e-09
50 8.5196068 -6.5897109 -6.4718776 985.38158 1.3996941e-08 7.0149406e-10 2.002272e-09
60 7.9273388 -6.5815192 -6.4718775 1000.4024 1.4000005e-08 3.5766629e-10 2.4944703e-09
70 7.2715879 -6.5724494 -6.4718773 1016.9932 1.3996503e-08 6.2731503e-10 1.7010533e-09
80 6.5722375 -6.5627766 -6.4718771 1034.6361 1.3973603e-08 3.1142917e-10 2.808524e-09
90 5.8505991 -6.5527956 -6.4718769 1052.7794 1.3983301e-08 3.9931135e-10 2.6118214e-09
100 5.128708 -6.542811 -6.4718767 1070.8561 1.395586e-08 2.3152413e-10 2.8742755e-09
110 4.4285344 -6.5331269 -6.4718766 1088.305 1.3938374e-08 4.2173005e-10 2.3059886e-09
120 3.7711361 -6.5240343 -6.4718764 1104.5919 1.3915264e-08 2.5458038e-10 1.4864012e-09
130 3.1757964 -6.5158002 -6.4718762 1119.2319 1.3858843e-08 5.7490448e-10 2.6191823e-09
140 2.6591997 -6.5086551 -6.4718761 1131.8095 1.3814891e-08 3.5434633e-10 2.2009364e-09
150 2.2347034 -6.5027839 -6.471876 1141.9961 1.3781115e-08 5.0639594e-10 2.9032558e-09
160 1.9117661 -6.4983173 -6.471876 1149.564 1.3734288e-08 3.1954962e-10 2.6097446e-09
170 1.6955808 -6.4953273 -6.471876 1154.3946 1.3682252e-08 3.5426781e-10 2.9605676e-09
180 1.586949 -6.4938249 -6.471876 1156.4812 1.363e-08 4.0804881e-10 2.1707904e-09
190 1.5824056 -6.4937621 -6.4718761 1155.925 1.3532637e-08 4.0767685e-10 3.0091462e-09
200 1.6745831 -6.4950371 -6.4718762 1152.926 1.3455927e-08 2.953369e-10 2.5029298e-09
210 1.8527803 -6.4975018 -6.4718763 1147.7684 1.335224e-08 3.5042319e-10 3.0550064e-09
220 2.1036825 -6.5009721 -6.4718764 1140.8026 1.3239176e-08 3.5988448e-10 2.6852683e-09
230 2.4121721 -6.5052389 -6.4718766 1132.4243 1.3090019e-08 3.5004036e-10 2.8838602e-09
240 2.7621668 -6.5100798 -6.4718767 1123.0538 1.2946525e-08 4.1216361e-10 2.1105916e-09
250 3.1374274 -6.5152701 -6.4718768 1113.1152 1.277789e-08 5.9848318e-10 2.3087106e-09
260 3.5222906 -6.5205932 -6.471877 1103.0171 1.2591089e-08 2.0080182e-10 1.6969069e-09
270 3.9022942 -6.5258491 -6.4718771 1093.1369 1.2432232e-08 4.2494727e-10 1.7375594e-09
280 4.2646753 -6.5308612 -6.4718772 1083.8072 1.2268238e-08 6.1239266e-10 1.7005135e-09
290 4.598736 -6.5354816 -6.4718772 1075.306 1.2181179e-08 4.9338341e-10 2.1326848e-09
300 4.896078 -6.5395941 -6.4718773 1067.85 1.2098274e-08 3.4564838e-10 2.4199891e-09
310 5.150715 -6.543116 -6.4718773 1061.5918 1.2184958e-08 4.2383299e-10 2.2243759e-09
320 5.3590742 -6.5459978 -6.4718773 1056.6189 1.2312948e-08 3.5194185e-10 1.3856935e-09
330 5.5199009 -6.5482222 -6.4718773 1052.9565 1.2573918e-08 4.2401322e-10 2.9882e-09
340 5.6340787 -6.5498013 -6.4718773 1050.5719 1.2821551e-08 5.8802825e-10 2.7333289e-09
350 5.7043792 -6.5507736 -6.4718772 1049.3813 1.3067314e-08 4.0014945e-10 2.3564728e-09
360 5.7351548 -6.5511992 -6.4718772 1049.2581 1.331283e-08 4.1684815e-10 1.735621e-09
370 5.7319891 -6.5511553 -6.4718771 1050.042 1.354018e-08 3.8495426e-10 2.4460056e-09
380 5.7013193 -6.5507311 -6.4718771 1051.5496 1.3734888e-08 3.5333605e-10 2.5174342e-09
390 5.6500487 -6.5500219 -6.471877 1053.5847 1.3892287e-08 3.8154957e-10 1.77358e-09
400 5.5851679 -6.5491245 -6.471877 1055.9489 1.3988171e-08 5.8769536e-10 1.9262201e-09
410 5.5134009 -6.5481319 -6.4718769 1058.4508 1.4088779e-08 3.6754739e-10 2.7586362e-09
420 5.4408957 -6.547129 -6.4718769 1060.9152 1.4139924e-08 4.9030281e-10 3.2871245e-09
430 5.3729707 -6.5461895 -6.4718768 1063.1898 1.4173041e-08 5.2345074e-10 3.5995984e-09
440 5.3139284 -6.5453729 -6.4718768 1065.1506 1.4191516e-08 5.9481094e-10 2.5005297e-09
450 5.2669383 -6.5447229 -6.4718768 1066.7054 1.4168424e-08 3.0799668e-10 2.0864191e-09
460 5.2339881 -6.5442672 -6.4718768 1067.7958 1.4163444e-08 6.3927736e-10 2.2872669e-09
470 5.2158979 -6.544017 -6.4718768 1068.3968 1.413819e-08 5.5108262e-10 4.4334972e-09
480 5.2123873 -6.5439685 -6.4718768 1068.5155 1.4083227e-08 3.9249548e-10 2.5568235e-09
490 5.2221849 -6.544104 -6.4718768 1068.188 1.4035287e-08 2.1988631e-10 2.1264034e-09
500 5.2431716 -6.5443943 -6.4718768 1067.4759 1.3968666e-08 3.9100701e-10 3.290368e-09
Loop time of 0.170182 on 4 procs for 500 steps with 108 atoms
Performance: 253.846 ns/day, 0.095 hours/ns, 2938.035 timesteps/s
99.7% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.0012069 | 0.0012994 | 0.0013512 | 0.2 | 0.76
Neigh | 0.0048233 | 0.0050244 | 0.0053881 | 0.3 | 2.95
Comm | 0.0072462 | 0.0078013 | 0.008175 | 0.4 | 4.58
Output | 0.0080632 | 0.0081244 | 0.0082899 | 0.1 | 4.77
Modify | 0.1476 | 0.14764 | 0.14768 | 0.0 | 86.75
Other | | 0.0002961 | | | 0.17
Nlocal: 27 ave 31 max 24 min
Histogram: 1 0 1 0 1 0 0 0 0 1
Nghost: 135 ave 138 max 131 min
Histogram: 1 0 0 0 0 1 0 1 0 1
Neighs: 162 ave 191 max 148 min
Histogram: 1 2 0 0 0 0 0 0 0 1
Total # of neighbors = 648
Ave neighs/atom = 6
Neighbor list builds = 500
Dangerous builds not checked
Total wall time: 0:00:00

175
examples/numdiff/log.28Jan2022.log.numdiff.g++.1
View File

@ -1,175 +0,0 @@
LAMMPS (7 Jan 2022)
# Numerical difference calculation
# of error in forces and virial stress
# adjustable parameters
variable nsteps index 500 # length of run
variable nthermo index 10 # thermo output interval
variable ndump index 500 # dump output interval
variable nlat index 3 # size of box
variable fdelta index 1.0e-4 # displacement size
variable vdelta index 1.0e-6 # strain size
variable temp index 10.0 # temperature
units metal
atom_style atomic
atom_modify map yes
lattice fcc 5.358000
Lattice spacing in x,y,z = 5.358 5.358 5.358
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
region box block 0 3 0 ${nlat} 0 ${nlat}
region box block 0 3 0 3 0 ${nlat}
region box block 0 3 0 3 0 3
create_box 1 box
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 108 atoms
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
create_atoms CPU = 0.000 seconds
mass 1 39.903
velocity all create ${temp} 2357 mom yes dist gaussian
velocity all create 10.0 2357 mom yes dist gaussian
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
timestep 0.001
fix nve all nve
# define numerical force calculation
fix numforce all numdiff ${nthermo} ${fdelta}
fix numforce all numdiff 10 ${fdelta}
fix numforce all numdiff 10 1.0e-4
variable ferrx atom f_numforce[1]-fx
variable ferry atom f_numforce[2]-fy
variable ferrz atom f_numforce[3]-fz
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
compute faverrsq all reduce ave v_ferrsq
variable fsq atom fx^2+fy^2+fz^2
compute favsq all reduce ave v_fsq
variable frelerr equal sqrt(c_faverrsq/c_favsq)
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
# define numerical virial stress tensor calculation
compute myvirial all pressure NULL virial
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
fix numvirial all numdiff/virial 10 ${vdelta}
fix numvirial all numdiff/virial 10 1.0e-6
variable errxx equal f_numvirial[1]-c_myvirial[1]
variable erryy equal f_numvirial[2]-c_myvirial[2]
variable errzz equal f_numvirial[3]-c_myvirial[3]
variable erryz equal f_numvirial[4]-c_myvirial[6]
variable errxz equal f_numvirial[5]-c_myvirial[5]
variable errxy equal f_numvirial[6]-c_myvirial[4]
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
variable vrelerr equal sqrt(v_verrsq/v_vsq)
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr
thermo ${nthermo}
thermo 10
run ${nsteps}
run 500
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.9407173
ghost atom cutoff = 5.9407173
binsize = 2.9703587, bins = 6 6 6
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.083 | 6.083 | 6.083 Mbytes
Step Temp PotEng TotEng Press v_frelerr v_vrelerr
0 10 -7.0259569 -6.8876486 28.564278 19203.344 1.5660292e-06
10 9.9376583 -7.0250947 -6.8876486 30.254762 1.5040965e-08 2.1991382e-07
20 9.7520139 -7.022527 -6.8876485 35.28505 1.4756358e-08 2.6265315e-06
30 9.4477557 -7.0183188 -6.8876485 43.519863 1.4688198e-08 2.6356166e-07
40 9.0330215 -7.0125826 -6.8876484 54.727797 1.4637921e-08 5.2292327e-08
50 8.5192918 -7.0054772 -6.8876483 68.585553 1.4587854e-08 7.1324716e-08
60 7.9212026 -6.997205 -6.8876481 84.684636 1.4525561e-08 3.1108149e-08
70 7.2562592 -6.9880081 -6.8876479 102.54088 1.450885e-08 3.2311094e-08
80 6.5444294 -6.9781627 -6.8876478 121.60715 1.4444738e-08 2.1776998e-08
90 5.8075961 -6.9679715 -6.8876476 141.2895 1.4493562e-08 2.0400898e-08
100 5.0688629 -6.957754 -6.8876474 160.9668 1.445455e-08 1.2636688e-08
110 4.3517145 -6.947835 -6.8876472 180.0135 1.4460371e-08 1.2528038e-08
120 3.6790589 -6.9385314 -6.887647 197.82486 1.4371757e-08 1.4489522e-08
130 3.0721984 -6.9301379 -6.8876468 213.84331 1.4364708e-08 1.2461922e-08
140 2.5497991 -6.9229125 -6.8876467 227.58429 1.4330926e-08 9.3913926e-09
150 2.1269443 -6.917064 -6.8876466 238.6596 1.4287002e-08 4.1510266e-09
160 1.8143642 -6.9127407 -6.8876465 246.79599 1.4282669e-08 7.7048281e-09
170 1.6179191 -6.9100237 -6.8876465 251.84748 1.42726e-08 1.2719973e-08
180 1.5383946 -6.9089239 -6.8876466 253.79991 1.4236534e-08 8.1200831e-09
190 1.5716287 -6.9093836 -6.8876467 252.76745 1.41706e-08 6.5670612e-09
200 1.7089493 -6.911283 -6.8876468 248.98142 1.4096463e-08 1.1685863e-08
210 1.9378716 -6.9144493 -6.8876469 242.77289 1.4008978e-08 1.1226902e-08
220 2.2429731 -6.9186692 -6.887647 234.55055 1.3886901e-08 9.9914102e-09
230 2.606862 -6.9237023 -6.8876472 224.77626 1.3864576e-08 1.1540228e-08
240 3.0111524 -6.9292941 -6.8876474 213.93996 1.3696314e-08 1.1697747e-08
250 3.4373794 -6.9351893 -6.8876475 202.53583 1.3626701e-08 1.0398197e-08
260 3.8678047 -6.9411426 -6.8876476 191.04084 1.3489489e-08 6.6603364e-09
270 4.2860853 -6.9469279 -6.8876478 179.89646 1.3312014e-08 1.1687917e-08
280 4.6777954 -6.9523457 -6.8876479 169.49404 1.3081144e-08 1.1336675e-08
290 5.030805 -6.9572282 -6.887648 160.16371 1.2947385e-08 1.7342825e-08
300 5.3355278 -6.9614428 -6.887648 152.16682 1.2893673e-08 1.7510534e-08
310 5.5850532 -6.964894 -6.887648 145.69148 1.2842022e-08 1.2782546e-08
320 5.7751794 -6.9675236 -6.8876481 140.85102 1.2903488e-08 1.5319437e-08
330 5.9043601 -6.9693103 -6.887648 137.68497 1.3076809e-08 1.1208999e-08
340 5.9735784 -6.9702676 -6.887648 136.16232 1.3296904e-08 1.891087e-08
350 5.9861549 -6.9704415 -6.887648 136.18679 1.3504051e-08 2.5783601e-08
360 5.947496 -6.9699067 -6.8876479 137.60397 1.3731112e-08 2.0556839e-08
370 5.8647874 -6.9687627 -6.8876478 140.2101 1.4009878e-08 2.1771736e-08
380 5.7466376 -6.9671285 -6.8876477 143.76234 1.4092054e-08 1.1085162e-08
390 5.6026773 -6.9651374 -6.8876477 147.99019 1.4282872e-08 2.0221602e-08
400 5.4431231 -6.9629305 -6.8876476 152.60787 1.4317739e-08 1.7076065e-08
410 5.2783192 -6.960651 -6.8876475 157.32722 1.4415075e-08 2.5031776e-08
420 5.1182723 -6.9584374 -6.8876474 161.87063 1.4441435e-08 2.2519289e-08
430 4.9722 -6.956417 -6.8876473 165.98344 1.4550624e-08 2.4512613e-08
440 4.8481153 -6.9547008 -6.8876473 169.44527 1.4544672e-08 1.4758301e-08
450 4.7524707 -6.9533779 -6.8876472 172.07964 1.4546492e-08 1.324687e-08
460 4.6898817 -6.9525122 -6.8876472 173.76132 1.4537475e-08 1.351367e-08
470 4.6629495 -6.9521397 -6.8876472 174.42109 1.4530458e-08 1.521106e-08
480 4.6721922 -6.9522675 -6.8876472 174.04742 1.4543785e-08 1.0905422e-08
490 4.7160887 -6.9528747 -6.8876473 172.68525 1.4545591e-08 2.0128525e-08
500 4.7912313 -6.953914 -6.8876473 170.43183 1.4438981e-08 1.6062775e-08
Loop time of 0.837333 on 1 procs for 500 steps with 108 atoms
Performance: 51.592 ns/day, 0.465 hours/ns, 597.134 timesteps/s
99.8% 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.0097726 | 0.0097726 | 0.0097726 | 0.0 | 1.17
Neigh | 0.03095 | 0.03095 | 0.03095 | 0.0 | 3.70
Comm | 0.005564 | 0.005564 | 0.005564 | 0.0 | 0.66
Output | 0.0042451 | 0.0042451 | 0.0042451 | 0.0 | 0.51
Modify | 0.78618 | 0.78618 | 0.78618 | 0.0 | 93.89
Other | | 0.0006258 | | | 0.07
Nlocal: 108 ave 108 max 108 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 558 ave 558 max 558 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 972 ave 972 max 972 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 972
Ave neighs/atom = 9
Neighbor list builds = 500
Dangerous builds not checked
Total wall time: 0:00:00

175
examples/numdiff/log.28Jan2022.log.numdiff.g++.4
View File

@ -1,175 +0,0 @@
LAMMPS (7 Jan 2022)
# Numerical difference calculation
# of error in forces and virial stress
# adjustable parameters
variable nsteps index 500 # length of run
variable nthermo index 10 # thermo output interval
variable ndump index 500 # dump output interval
variable nlat index 3 # size of box
variable fdelta index 1.0e-4 # displacement size
variable vdelta index 1.0e-6 # strain size
variable temp index 10.0 # temperature
units metal
atom_style atomic
atom_modify map yes
lattice fcc 5.358000
Lattice spacing in x,y,z = 5.358 5.358 5.358
region box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
region box block 0 3 0 ${nlat} 0 ${nlat}
region box block 0 3 0 3 0 ${nlat}
region box block 0 3 0 3 0 3
create_box 1 box
Created orthogonal box = (0 0 0) to (16.074 16.074 16.074)
1 by 2 by 2 MPI processor grid
create_atoms 1 box
Created 108 atoms
using lattice units in orthogonal box = (0 0 0) to (16.074 16.074 16.074)
create_atoms CPU = 0.000 seconds
mass 1 39.903
velocity all create ${temp} 2357 mom yes dist gaussian
velocity all create 10.0 2357 mom yes dist gaussian
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42
neighbor 0.0 bin
neigh_modify every 1 delay 0 check no
timestep 0.001
fix nve all nve
# define numerical force calculation
fix numforce all numdiff ${nthermo} ${fdelta}
fix numforce all numdiff 10 ${fdelta}
fix numforce all numdiff 10 1.0e-4
variable ferrx atom f_numforce[1]-fx
variable ferry atom f_numforce[2]-fy
variable ferrz atom f_numforce[3]-fz
variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
compute faverrsq all reduce ave v_ferrsq
variable fsq atom fx^2+fy^2+fz^2
compute favsq all reduce ave v_fsq
variable frelerr equal sqrt(c_faverrsq/c_favsq)
dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz
dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz
# define numerical virial stress tensor calculation
compute myvirial all pressure NULL virial
fix numvirial all numdiff/virial ${nthermo} ${vdelta}
fix numvirial all numdiff/virial 10 ${vdelta}
fix numvirial all numdiff/virial 10 1.0e-6
variable errxx equal f_numvirial[1]-c_myvirial[1]
variable erryy equal f_numvirial[2]-c_myvirial[2]
variable errzz equal f_numvirial[3]-c_myvirial[3]
variable erryz equal f_numvirial[4]-c_myvirial[6]
variable errxz equal f_numvirial[5]-c_myvirial[5]
variable errxy equal f_numvirial[6]-c_myvirial[4]
variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2"
variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2"
variable vrelerr equal sqrt(v_verrsq/v_vsq)
thermo_style custom step temp pe etotal press v_frelerr v_vrelerr
thermo ${nthermo}
thermo 10
run ${nsteps}
run 500
generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 0 steps, check no
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.9407173
ghost atom cutoff = 5.9407173
binsize = 2.9703587, bins = 6 6 6
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cubic, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.067 | 6.067 | 6.067 Mbytes
Step Temp PotEng TotEng Press v_frelerr v_vrelerr
0 10 -7.0259569 -6.8876486 28.564278 10664.391 9.1481187e-08
10 9.9388179 -7.0251107 -6.8876486 30.21736 1.4771865e-08 1.3452512e-07
20 9.7572185 -7.022599 -6.8876485 35.123527 1.437525e-08 8.0966999e-07
30 9.4606673 -7.0184974 -6.8876484 43.132052 1.4375468e-08 1.990012e-08
40 9.0579092 -7.0129268 -6.8876483 54.000927 1.4350331e-08 1.7239028e-08
50 8.5607685 -7.0060508 -6.8876482 67.403151 1.4353284e-08 7.813181e-09
60 7.9838726 -6.9980717 -6.8876481 82.935358 1.4398078e-08 2.022316e-08
70 7.3442875 -6.9892255 -6.8876479 100.12892 1.434409e-08 7.5938179e-09
80 6.6610579 -6.9797757 -6.8876477 118.46358 1.4324787e-08 7.1972571e-09
90 5.9546462 -6.9700053 -6.8876476 137.38365 1.4322718e-08 4.3978378e-09
100 5.2462727 -6.9602077 -6.8876474 156.31651 1.4273172e-08 4.6728038e-09
110 4.5571706 -6.9506767 -6.8876472 174.69294 1.4266163e-08 3.522225e-09
120 3.9077807 -6.9416949 -6.887647 191.96859 1.42241e-08 3.5357511e-09
130 3.3169241 -6.9335227 -6.8876469 207.64566 1.4203813e-08 2.5182488e-09
140 2.8010028 -6.926387 -6.8876468 221.29333 1.4164215e-08 2.3112509e-09
150 2.3732854 -6.9204712 -6.8876467 232.5658 1.4134122e-08 1.9368963e-09
160 2.0433329 -6.9159076 -6.8876466 241.21647 1.4095473e-08 3.6806452e-09
170 1.8166184 -6.912772 -6.8876466 247.10715 1.4049531e-08 2.5319322e-09
180 1.6943727 -6.9110813 -6.8876467 250.21143 1.3997912e-08 1.952404e-09
190 1.6736731 -6.910795 -6.8876467 250.61203 1.3915487e-08 1.4271767e-09
200 1.7477659 -6.9118199 -6.8876468 248.49223 1.3850618e-08 2.4515718e-09
210 1.9065921 -6.9140167 -6.8876469 244.12226 1.3747916e-08 1.7957531e-09
220 2.1374676 -6.91721 -6.887647 237.84173 1.3674779e-08 2.6613511e-09
230 2.4258607 -6.9211989 -6.8876472 230.0395 1.3565712e-08 3.0777067e-09
240 2.7562034 -6.9257679 -6.8876473 221.13265 1.3440442e-08 1.7111501e-09
250 3.1126827 -6.9306984 -6.8876474 211.54566 1.3270914e-08 1.6690112e-09
260 3.4799641 -6.9357784 -6.8876476 201.69126 1.3105092e-08 2.1637558e-09
270 3.8438148 -6.9408108 -6.8876477 191.95361 1.2962846e-08 4.4613506e-09
280 4.191607 -6.9456212 -6.8876478 182.6745 1.2846383e-08 3.3730203e-09
290 4.5126922 -6.9500621 -6.8876478 174.14285 1.2710692e-08 2.2809889e-09
300 4.7986487 -6.9540172 -6.8876479 166.58747 1.2657778e-08 6.9880891e-09
310 5.0434083 -6.9574025 -6.8876479 160.17316 1.266381e-08 4.2486217e-09
320 5.243275 -6.9601668 -6.8876479 154.99974 1.279856e-08 5.1505673e-09
330 5.3968455 -6.9622908 -6.8876479 151.1038 1.2981831e-08 3.3339333e-09
340 5.5048468 -6.9637845 -6.8876479 148.46296 1.3159021e-08 1.7881393e-09
350 5.569902 -6.9646843 -6.8876479 147.00205 1.3439465e-08 5.6876219e-09
360 5.5962378 -6.9650485 -6.8876478 146.60113 1.3645943e-08 7.233847e-09
370 5.5893468 -6.9649531 -6.8876478 147.10471 1.3829581e-08 4.5514318e-09
380 5.5556199 -6.9644866 -6.8876477 148.33195 1.4005893e-08 4.2971846e-09
390 5.5019639 -6.9637444 -6.8876476 150.08725 1.4157454e-08 3.3564262e-09
400 5.4354239 -6.962824 -6.8876476 152.17073 1.4226422e-08 4.2393923e-09
410 5.3628267 -6.9618199 -6.8876475 154.38825 1.4302679e-08 3.8937698e-09
420 5.2904639 -6.960819 -6.8876475 156.56034 1.4381099e-08 4.315875e-09
430 5.2238282 -6.9598973 -6.8876474 158.52969 1.4426567e-08 4.2658185e-09
440 5.1674149 -6.9591171 -6.8876474 160.16704 1.4453381e-08 5.7055268e-09
450 5.1245913 -6.9585248 -6.8876474 161.37513 1.4449488e-08 4.4308801e-09
460 5.0975361 -6.9581506 -6.8876474 162.09077 1.4445596e-08 5.8269923e-09
470 5.0872416 -6.9580082 -6.8876474 162.28517 1.4444348e-08 4.8263194e-09
480 5.0935712 -6.9580957 -6.8876474 161.96268 1.4411666e-08 6.222228e-09
490 5.115362 -6.9583971 -6.8876474 161.15816 1.4369716e-08 3.3926077e-09
500 5.1505605 -6.958884 -6.8876474 159.9333 1.4288515e-08 3.8845251e-09
Loop time of 0.252598 on 4 procs for 500 steps with 108 atoms
Performance: 171.023 ns/day, 0.140 hours/ns, 1979.430 timesteps/s
99.8% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.0021545 | 0.0022845 | 0.0023794 | 0.2 | 0.90
Neigh | 0.0063887 | 0.0067604 | 0.0069752 | 0.3 | 2.68
Comm | 0.01048 | 0.010916 | 0.011408 | 0.3 | 4.32
Output | 0.0026603 | 0.0027399 | 0.0029738 | 0.3 | 1.08
Modify | 0.2295 | 0.22952 | 0.22954 | 0.0 | 90.86
Other | | 0.0003814 | | | 0.15
Nlocal: 27 ave 29 max 25 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Nghost: 325 ave 327 max 323 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Neighs: 243 ave 273 max 228 min
Histogram: 1 1 1 0 0 0 0 0 0 1
Total # of neighbors = 972
Ave neighs/atom = 9
Neighbor list builds = 500
Dangerous builds not checked
Total wall time: 0:00:00

49
lib/atc/LammpsInterface.cpp
View File

@ -387,7 +387,7 @@ double LammpsInterface::atom_quantity_conversion(FundamentalAtomQuantity quantit
int LammpsInterface::dimension() const { return lammps_->domain->dimension; }
int LammpsInterface::nregion() const { return lammps_->domain->nregion; }
int LammpsInterface::nregion() const { return lammps_->domain->get_region_list().size(); }
void LammpsInterface::box_bounds(double & boxxlo, double & boxxhi,
double & boxylo, double & boxyhi,
@ -527,14 +527,15 @@ void LammpsInterface::box_periodicity(int & xperiodic,
zperiodic = lammps_->domain->zperiodic;
}
int LammpsInterface::region_id(const char * regionName) const {
int nregion = this->nregion();
for (int iregion = 0; iregion < nregion; iregion++) {
if (strcmp(regionName, region_name(iregion)) == 0) {
int LammpsInterface::region_id(const char *regionName) const {
auto regions = lammps_->domain->get_region_list();
int iregion = 0;
for (auto reg : regions) {
if (strcmp(regionName, reg->id) == 0) {
return iregion;
}
++iregion;
}
throw ATC_Error("Region has not been defined");
return -1;
}
@ -1322,61 +1323,73 @@ int** LammpsInterface::bond_list() const { return lammps_->neighbor->bondlist;
char * LammpsInterface::region_name(int iRegion) const
{
return lammps_->domain->regions[iRegion]->id;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->id;
}
char * LammpsInterface::region_style(int iRegion) const
{
return lammps_->domain->regions[iRegion]->style;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->style;
}
double LammpsInterface::region_xlo(int iRegion) const
{
return lammps_->domain->regions[iRegion]->extent_xlo;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->extent_xlo;
}
double LammpsInterface::region_xhi(int iRegion) const
{
return lammps_->domain->regions[iRegion]->extent_xhi;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->extent_xhi;
}
double LammpsInterface::region_ylo(int iRegion) const
{
return lammps_->domain->regions[iRegion]->extent_ylo;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->extent_ylo;
}
double LammpsInterface::region_yhi(int iRegion) const
{
return lammps_->domain->regions[iRegion]->extent_yhi;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->extent_yhi;
}
double LammpsInterface::region_zlo(int iRegion) const
{
return lammps_->domain->regions[iRegion]->extent_zlo;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->extent_zlo;
}
double LammpsInterface::region_zhi(int iRegion) const
{
return lammps_->domain->regions[iRegion]->extent_zhi;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->extent_zhi;
}
double LammpsInterface::region_xscale(int iRegion) const
{
return lammps_->domain->regions[iRegion]->xscale;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->xscale;
}
double LammpsInterface::region_yscale(int iRegion) const
{
return lammps_->domain->regions[iRegion]->yscale;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->yscale;
}
double LammpsInterface::region_zscale(int iRegion) const
{
return lammps_->domain->regions[iRegion]->zscale;
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->zscale;
}
int LammpsInterface::region_match(int iRegion, double x, double y, double z) const {
return lammps_->domain->regions[iRegion]->match(x,y,z);
auto regions = lammps_->domain->get_region_list();
return regions[iRegion]->match(x,y,z);
}
// -----------------------------------------------------------------

4
src/.gitignore vendored
View File

@ -453,6 +453,8 @@
/compute_basal_atom.h
/compute_body_local.cpp
/compute_body_local.h
/compute_born_matrix.cpp
/compute_born_matrix.h
/compute_cnp_atom.cpp
/compute_cnp_atom.h
/compute_damage_atom.cpp
@ -1091,6 +1093,8 @@
/pair_hdnnp.h
/pair_ilp_graphene_hbn.cpp
/pair_ilp_graphene_hbn.h
/pair_ilp_graphene_hbn_opt.cpp
/pair_ilp_graphene_hbn_opt.h
/pair_ilp_tmd.cpp
/pair_ilp_tmd.h
/pair_kolmogorov_crespi_full.cpp

12
src/ADIOS/dump_atom_adios.cpp
View File

@ -134,10 +134,10 @@ void DumpAtomADIOS::write()
// Now we know the global size and the local subset size and offset
// of the atoms table
auto nAtomsGlobal = static_cast<size_t>(ntotal);
auto startRow = static_cast<size_t>(atomOffset);
auto nAtomsLocal = static_cast<size_t>(nme);
auto nColumns = static_cast<size_t>(size_one);
auto nAtomsGlobal = static_cast<size_t>(ntotal);
auto startRow = static_cast<size_t>(atomOffset);
auto nAtomsLocal = static_cast<size_t>(nme);
auto nColumns = static_cast<size_t>(size_one);
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
@ -238,7 +238,7 @@ void DumpAtomADIOS::init_style()
columnNames = {"id", "type", "xs", "ys", "zs", "ix", "iy", "iz"};
}
for (int icol = 0; icol < (int)columnNames.size(); ++icol)
for (int icol = 0; icol < (int) columnNames.size(); ++icol)
if (keyword_user[icol].size()) columnNames[icol] = keyword_user[icol];
// setup function ptrs
@ -296,7 +296,7 @@ void DumpAtomADIOS::init_style()
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
auto nColumns = static_cast<size_t>(size_one);
auto nColumns = static_cast<size_t>(size_one);
internal->io.DefineAttribute<std::string>("columns", columnNames.data(), nColumns);
internal->io.DefineAttribute<std::string>("columnstr", columns);
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);

14
src/ADIOS/dump_atom_adios.h
View File

@ -44,17 +44,3 @@ class DumpAtomADIOS : public DumpAtom {
#endif
#endif
/* ERROR/WARNING messages:
E: Cannot open dump file %s
The output file for the dump command cannot be opened. Check that the
path and name are correct.
E: Too much per-proc info for dump
Number of local atoms times number of columns must fit in a 32-bit
integer for dump.
*/

42
src/ADIOS/dump_custom_adios.cpp
View File

@ -146,10 +146,10 @@ void DumpCustomADIOS::write()
// Now we know the global size and the local subset size and offset
// of the atoms table
auto nAtomsGlobal = static_cast<size_t>(ntotal);
auto startRow = static_cast<size_t>(atomOffset);
auto nAtomsLocal = static_cast<size_t>(nme);
auto nColumns = static_cast<size_t>(size_one);
auto nAtomsGlobal = static_cast<size_t>(ntotal);
auto startRow = static_cast<size_t>(atomOffset);
auto nAtomsLocal = static_cast<size_t>(nme);
auto nColumns = static_cast<size_t>(size_one);
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
@ -221,8 +221,10 @@ void DumpCustomADIOS::init_style()
int icol = 0;
for (auto item : utils::split_words(columns_default)) {
if (combined.size()) combined += " ";
if (keyword_user[icol].size()) combined += keyword_user[icol];
else combined += item;
if (keyword_user[icol].size())
combined += keyword_user[icol];
else
combined += item;
++icol;
}
columns = utils::strdup(combined);
@ -249,34 +251,30 @@ void DumpCustomADIOS::init_style()
*/
// find current ptr for each compute,fix,variable
// check that fix frequency is acceptable
int icompute;
for (int i = 0; i < ncompute; i++) {
icompute = modify->find_compute(id_compute[i]);
if (icompute < 0) error->all(FLERR, "Could not find dump custom compute ID");
compute[i] = modify->compute[icompute];
compute[i] = modify->get_compute_by_id(id_compute[i]);
if (!compute[i])
error->all(FLERR, "Could not find dump custom/adios compute ID {}", id_compute[i]);
}
int ifix;
for (int i = 0; i < nfix; i++) {
ifix = modify->find_fix(id_fix[i]);
if (ifix < 0) error->all(FLERR, "Could not find dump custom fix ID");
fix[i] = modify->fix[ifix];
if (nevery % modify->fix[ifix]->peratom_freq)
error->all(FLERR, "Dump custom and fix not computed at compatible times");
fix[i] = modify->get_fix_by_id(id_fix[i]);
if (!fix[i]) error->all(FLERR, "Could not find dump custom/adios fix ID {}", id_fix[i]);
if (nevery % fix[i]->peratom_freq)
error->all(FLERR, "dump custom/adios and fix {} with ID {} not computed at compatible times",
fix[i]->style, id_fix[i]);
}
int ivariable;
for (int i = 0; i < nvariable; i++) {
ivariable = input->variable->find(id_variable[i]);
if (ivariable < 0) error->all(FLERR, "Could not find dump custom variable name");
if (ivariable < 0) error->all(FLERR, "Could not find dump custom/adios variable name");
variable[i] = ivariable;
}
// set index and check validity of region
if (iregion >= 0) {
iregion = domain->find_region(idregion);
if (iregion == -1) error->all(FLERR, "Region ID for dump custom does not exist");
}
if (idregion && !domain->get_region_by_id(idregion))
error->all(FLERR, "Region {} for dump custom/adios does not exist", idregion);
/* Define the group of variables for the atom style here since it's a fixed
* set */
@ -316,7 +314,7 @@ void DumpCustomADIOS::init_style()
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
auto nColumns = static_cast<size_t>(size_one);
auto nColumns = static_cast<size_t>(size_one);
internal->io.DefineAttribute<std::string>("columns", internal->columnNames.data(), nColumns);
internal->io.DefineAttribute<std::string>("columnstr", columns);
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);

40
src/ADIOS/dump_custom_adios.h
View File

@ -43,43 +43,3 @@ class DumpCustomADIOS : public DumpCustom {
#endif
#endif
/* ERROR/WARNING messages:
E: Cannot open dump file %s
The output file for the dump command cannot be opened. Check that the
path and name are correct.
E: Too much per-proc info for dump
Number of local atoms times number of columns must fit in a 32-bit
integer for dump.
E: Dump_modify format string is too short
There are more fields to be dumped in a line of output than your
format string specifies.
E: Could not find dump custom compute ID
Self-explanatory.
E: Could not find dump custom fix ID
Self-explanatory.
E: Dump custom and fix not computed at compatible times
The fix must produce per-atom quantities on timesteps that dump custom
needs them.
E: Could not find dump custom variable name
Self-explanatory.
E: Region ID for dump custom does not exist
Self-explanatory.
*/

4
src/ADIOS/reader_adios.cpp
View File

@ -221,7 +221,7 @@ bigint ReaderADIOS::read_header(double box[3][3], int &boxinfo, int &triclinic,
uint64_t rem = nAtomsTotal % comm->nprocs;
nAtoms = nAtomsTotal / comm->nprocs;
atomOffset = comm->me * nAtoms;
if (comm->me < (int)rem) {
if (comm->me < (int) rem) {
++nAtoms;
atomOffset += comm->me;
} else {
@ -421,7 +421,7 @@ void ReaderADIOS::read_atoms(int n, int nfield, double **fields)
adios2::Variable<double> varAtoms = internal->io.InquireVariable<double>("atoms");
if ((uint64_t)n != nAtoms)
if ((uint64_t) n != nAtoms)
error->one(FLERR,
"ReaderADIOS::read_atoms() expects 'n={}' equal to the number of "
"atoms (={}) for process {} in ADIOS file {}.",

4
src/ADIOS/reader_adios.h
View File

@ -40,8 +40,8 @@ class ReaderADIOS : public Reader {
int read_time(bigint &) override;
void skip() override;
bigint read_header(double[3][3], int &, int &, int, int, int *, char **, int, int, int &,
int &, int &, int &) override;
bigint read_header(double[3][3], int &, int &, int, int, int *, char **, int, int, int &, int &,
int &, int &) override;
void read_atoms(int, int, double **) override;
void open_file(const std::string &) override;

18
src/ASPHERE/compute_erotate_asphere.h
View File

@ -41,21 +41,3 @@ class ComputeERotateAsphere : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute erotate/asphere requires atom style ellipsoid or line or tri
Self-explanatory.
E: Compute erotate/asphere requires extended particles
This compute cannot be used with point particles.
*/

40
src/ASPHERE/compute_temp_asphere.h
View File

@ -52,43 +52,3 @@ class ComputeTempAsphere : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute temp/asphere requires atom style ellipsoid
Self-explanatory.
E: Compute temp/asphere requires extended particles
This compute cannot be used with point particles.
E: Could not find compute ID for temperature bias
Self-explanatory.
E: Bias compute does not calculate temperature
The specified compute must compute temperature.
E: Bias compute does not calculate a velocity bias
The specified compute must compute a bias for temperature.
E: Bias compute group does not match compute group
The specified compute must operate on the same group as the parent
compute.
E: Temperature compute degrees of freedom < 0
This should not happen if you are calculating the temperature
on a valid set of atoms.
*/

13
src/ASPHERE/fix_nh_asphere.h
View File

@ -35,16 +35,3 @@ class FixNHAsphere : public FixNH {
} // namespace LAMMPS_NS
#endif
/* ERROR/WARNING messages:
E: Compute nvt/nph/npt asphere requires atom style ellipsoid
Self-explanatory.
E: Fix nvt/nph/npt asphere requires extended particles
The shape setting for a particle in the fix group has shape = 0.0,
which means it is a point particle.
*/

12
src/ASPHERE/fix_nph_asphere.h
View File

@ -33,15 +33,3 @@ class FixNPHAsphere : public FixNHAsphere {
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control can not be used with fix nph/asphere
Self-explanatory.
E: Pressure control must be used with fix nph/asphere
Self-explanatory.
*/

12
src/ASPHERE/fix_npt_asphere.h
View File

@ -33,15 +33,3 @@ class FixNPTAsphere : public FixNHAsphere {
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control must be used with fix npt/asphere
Self-explanatory.
E: Pressure control must be used with fix npt/asphere
Self-explanatory.
*/

12
src/ASPHERE/fix_nve_asphere.h
View File

@ -39,15 +39,3 @@ class FixNVEAsphere : public FixNVE {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Compute nve/asphere requires atom style ellipsoid
Self-explanatory.
E: Fix nve/asphere requires extended particles
This fix can only be used for particles with a shape setting.
*/

18
src/ASPHERE/fix_nve_asphere_noforce.h
View File

@ -39,21 +39,3 @@ class FixNVEAsphereNoforce : public FixNVENoforce {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix nve/asphere/noforce requires atom style ellipsoid
Self-explanatory.
E: Fix nve/asphere/noforce requires extended particles
One of the particles is not an ellipsoid.
*/

22
src/ASPHERE/fix_nve_line.h
View File

@ -41,25 +41,3 @@ class FixNVELine : public FixNVE {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix nve/line requires atom style line
Self-explanatory.
E: Fix nve/line can only be used for 2d simulations
Self-explanatory.
E: Fix nve/line requires line particles
Self-explanatory.
*/

22
src/ASPHERE/fix_nve_tri.h
View File

@ -41,25 +41,3 @@ class FixNVETri : public FixNVE {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix nve/tri requires atom style tri
Self-explanatory.
E: Fix nve/tri can only be used for 3d simulations
Self-explanatory.
E: Fix nve/tri requires tri particles
Self-explanatory.
*/

12
src/ASPHERE/fix_nvt_asphere.h
View File

@ -33,15 +33,3 @@ class FixNVTAsphere : public FixNHAsphere {
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control must be used with fix nvt/asphere
Self-explanatory.
E: Pressure control can not be used with fix nvt/asphere
Self-explanatory.
*/

31
src/ASPHERE/pair_gayberne.h
View File

@ -72,34 +72,3 @@ class PairGayBerne : public Pair {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair gayberne requires atom style ellipsoid
Self-explanatory.
E: Pair gayberne requires atoms with same type have same shape
Self-explanatory.
E: Pair gayberne epsilon a,b,c coeffs are not all set
Each atom type involved in pair_style gayberne must
have these 3 coefficients set at least once.
E: Bad matrix inversion in mldivide3
This error should not occur unless the matrix is badly formed.
*/

23
src/ASPHERE/pair_line_lj.h
View File

@ -62,26 +62,3 @@ class PairLineLJ : public Pair {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair line/lj requires atom style line
Self-explanatory.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/

4
src/ASPHERE/pair_resquared.cpp
View File

@ -317,10 +317,10 @@ void PairRESquared::coeff(int narg, char **arg)
void PairRESquared::init_style()
{
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
if (!avec) error->all(FLERR, "Pair resquared requires atom style ellipsoid");
neighbor->add_request(this,NeighConst::REQ_DEFAULT);
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
// per-type shape precalculations
// require that atom shapes are identical within each type

34
src/ASPHERE/pair_resquared.h
View File

@ -92,37 +92,3 @@ class PairRESquared : public Pair {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair resquared requires atom style ellipsoid
Self-explanatory.
E: Pair resquared requires atoms with same type have same shape
Self-explanatory.
E: Pair resquared epsilon a,b,c coeffs are not all set
Self-explanatory.
E: Pair resquared epsilon and sigma coeffs are not all set
Self-explanatory.
E: Bad matrix inversion in mldivide3
This error should not occur unless the matrix is badly formed.
*/

18
src/ASPHERE/pair_tri_lj.h
View File

@ -60,21 +60,3 @@ class PairTriLJ : public Pair {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair tri/lj requires atom style tri
Self-explanatory.
*/

78
src/AWPMD/atom_vec_wavepacket.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -20,8 +19,6 @@
#include "atom.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -32,33 +29,29 @@ AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp)
molecular = Atom::ATOMIC;
forceclearflag = 1;
atom->wavepacket_flag = 1;
atom->electron_flag = 1; // compatible with eff
atom->q_flag = atom->spin_flag = atom->eradius_flag =
atom->ervel_flag = atom->erforce_flag = 1;
atom->cs_flag = atom->csforce_flag =
atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
atom->wavepacket_flag = atom->q_flag = atom->spin_flag = atom->eradius_flag = 1;
atom->ervel_flag = atom->erforce_flag = atom->cs_flag = atom->csforce_flag = 1;
atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = (char *)
"q spin eradius ervel erforce cs csforce "
"vforce ervelforce etag";
fields_copy = (char *) "q spin eradius ervel cs etag";
fields_comm = (char *) "eradius";
fields_comm_vel = (char *) "eradius ervel cs";
fields_reverse = (char *) "erforce ervelforce vforce csforce";
fields_border = (char *) "q spin eradius etag";
fields_border_vel = (char *) "q spin eradius etag ervel cs";
fields_exchange = (char *) "q spin eradius ervel etag cs";
fields_restart = (char *) "q spin eradius ervel etag cs";
fields_create = (char *) "q spin eradius ervel etag cs";
fields_data_atom = (char *) "id type q spin eradius etag cs x";
fields_data_vel = (char *) "id v ervel";
fields_grow = {"q", "spin", "eradius", "ervel", "erforce",
"cs", "csforce", "vforce", "ervelforce", "etag"};
fields_copy = {"q", "spin", "eradius", "ervel", "cs", "etag"};
fields_comm = {"eradius"};
fields_comm_vel = {"eradius", "ervel", "cs"};
fields_reverse = {"erforce", "ervelforce", "vforce", "csforce"};
fields_border = {"q", "spin", "eradius", "etag"};
fields_border_vel = {"q", "spin", "eradius", "etag", "ervel", "cs"};
fields_exchange = {"q", "spin", "eradius", "ervel", "etag", "cs"};
fields_restart = {"q", "spin", "eradius", "ervel", "etag", "cs"};
fields_create = {"q", "spin", "eradius", "ervel", "etag", "cs"};
fields_data_atom = {"id", "type", "q", "spin", "eradius", "etag", "cs", "x"};
fields_data_vel = {"id", "v", "ervel"};
setup_fields();
}
@ -84,7 +77,7 @@ void AtomVecWavepacket::grow_pointers()
void AtomVecWavepacket::force_clear(int n, size_t nbytes)
{
memset(&erforce[n],0,nbytes);
memset(&erforce[n], 0, nbytes);
}
/* ----------------------------------------------------------------------
@ -112,12 +105,12 @@ void AtomVecWavepacket::data_atom_post(int ilocal)
return -1 if name is unknown to this atom style
------------------------------------------------------------------------- */
int AtomVecWavepacket::property_atom(char *name)
int AtomVecWavepacket::property_atom(const std::string &name)
{
if (strcmp(name,"spin") == 0) return 0;
if (strcmp(name,"eradius") == 0) return 1;
if (strcmp(name,"ervel") == 0) return 2;
if (strcmp(name,"erforce") == 0) return 3;
if (name == "spin") return 0;
if (name == "eradius") return 1;
if (name == "ervel") return 2;
if (name == "erforce") return 3;
return -1;
}
@ -126,34 +119,41 @@ int AtomVecWavepacket::property_atom(char *name)
index maps to data specific to this atom style
------------------------------------------------------------------------- */
void AtomVecWavepacket::pack_property_atom(int index, double *buf,
int nvalues, int groupbit)
void AtomVecWavepacket::pack_property_atom(int index, double *buf, int nvalues, int groupbit)
{
int nlocal = atom->nlocal;
int n = 0;
if (index == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = spin[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = spin[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = eradius[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = eradius[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = ervel[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = ervel[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 3) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = erforce[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = erforce[i];
else
buf[n] = 0.0;
n += nvalues;
}
}

2
src/AWPMD/atom_vec_wavepacket.h
View File

@ -32,7 +32,7 @@ class AtomVecWavepacket : public AtomVec {
void force_clear(int, size_t) override;
void create_atom_post(int) override;
void data_atom_post(int) override;
int property_atom(char *) override;
int property_atom(const std::string &) override;
void pack_property_atom(int, double *, int, int) override;
private:

39
src/BOCS/compute_pressure_bocs.h
View File

@ -80,42 +80,3 @@ class ComputePressureBocs : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute pressure must use group all
Virial contributions computed by potentials (pair, bond, etc) are
computed on all atoms.
E: Could not find compute pressure temperature ID
The compute ID for calculating temperature does not exist.
E: Compute pressure temperature ID does not compute temperature
The compute ID assigned to a pressure computation must compute
temperature.
E: Compute pressure requires temperature ID to include kinetic energy
The keflag cannot be used unless a temperature compute is provided.
E: Virial was not tallied on needed timestep
You are using a thermo keyword that requires potentials to
have tallied the virial, but they didn't on this timestep. See the
variable doc page for ideas on how to make this work.
E: Must use 'kspace_modify pressure/scalar no' for tensor components with kspace_style msm
Otherwise MSM will compute only a scalar pressure. See the kspace_modify
command for details on this setting.
*/

135
src/BOCS/fix_bocs.h
View File

@ -30,7 +30,7 @@ namespace LAMMPS_NS {
class FixBocs : public Fix {
public:
FixBocs(class LAMMPS *, int, char **); // MRD NJD
~FixBocs() override; // MRD NJD
~FixBocs() override; // MRD NJD
int setmask() override;
void init() override;
void setup(int) override;
@ -169,136 +169,3 @@ class FixBocs : public Fix {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: CG basis type XXX is not recognized
See second line of message for supported basis types.
E: Target temperature for fix bocs cannot be 0.0
Self-explanatory.
E: Invalid fix bocs command for a 2d simulation
Cannot control z dimension in a 2d model.
E: Fix bocs dilate group ID does not exist
Self-explanatory.
E: Invalid fix bocs command pressure settings
If multiple dimensions are coupled, those dimensions must be
specified.
E: Cannot use fix bocs on a non-periodic dimension
When specifying a diagonal pressure component, the dimension must be
periodic.
E: Cannot use fix bocs on a 2nd non-periodic dimension
When specifying an off-diagonal pressure component, the 2nd of the two
dimensions must be periodic. E.g. if the xy component is specified,
then the y dimension must be periodic.
E: Cannot use fix bocs with yz scaling when z is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with xz scaling when z is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with xy scaling when y is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with both yz dynamics and yz scaling
Self-explanatory.
E: Cannot use fix bocs with both xz dynamics and xz scaling
Self-explanatory.
E: Cannot use fix bocs with both xy dynamics and xy scaling
Self-explanatory.
E: Can not specify Pxy/Pxz/Pyz in fix bocs with non-triclinic box
Only triclinic boxes can be used with off-diagonal pressure components.
See the region prism command for details.
E: Invalid fix bocs pressure settings
Settings for coupled dimensions must be the same.
E: Using update dipole flag requires atom style sphere
Self-explanatory.
E: Using update dipole flag requires atom attribute mu
Self-explanatory.
E: The dlm flag must be used with update dipole
Self-explanatory.
E: Fix bocs damping parameters must be > 0.0
Self-explanatory.
E: Cannot use fix npt and fix deform on same component of stress tensor
This would be changing the same box dimension twice.
E: Temperature ID for fix bocs does not exist
Self-explanatory.
E: Pressure ID for fix bocs does not exist
Self-explanatory.
E: Fix bocs has tilted box too far in one step - periodic cell is too far from equilibrium state
Self-explanatory. The change in the box tilt is too extreme
on a short timescale.
E: Could not find fix_modify temperature ID
The compute ID for computing temperature does not exist.
E: Fix_modify temperature ID does not compute temperature
The compute ID assigned to the fix must compute temperature.
W: Temperature for fix modify is not for group all
The temperature compute is being used with a pressure calculation
which does operate on group all, so this may be inconsistent.
E: Pressure ID for fix modify does not exist
Self-explanatory.
E: Could not find fix_modify pressure ID
The compute ID for computing pressure does not exist.
E: Fix_modify pressure ID does not compute pressure
The compute ID assigned to the fix must compute pressure.
*/

24
src/BODY/body_nparticle.h
View File

@ -53,27 +53,3 @@ class BodyNparticle : public Body {
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body nparticle command
Arguments in atom-style command are not correct.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body nparticle
Eigensolve for rigid body was not sufficiently accurate.
*/

29
src/BODY/body_rounded_polygon.h
View File

@ -57,32 +57,3 @@ class BodyRoundedPolygon : public Body {
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body rounded/polygon command
Arguments in atom-style command are not correct.
E: Invalid format in Bodies section of data file
The specified number of integer or floating point values does not
appear.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body nparticle
Eigensolve for rigid body was not sufficiently accurate.
*/

29
src/BODY/body_rounded_polyhedron.h
View File

@ -59,32 +59,3 @@ class BodyRoundedPolyhedron : public Body {
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body rounded/polyhedron command
Arguments in atom-style command are not correct.
E: Invalid format in Bodies section of data file
The specified number of integer or floating point values does not
appear.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body rounded/polyhedron
Eigensolve for rigid body was not sufficiently accurate.
*/

22
src/BODY/compute_body_local.h
View File

@ -49,25 +49,3 @@ class ComputeBodyLocal : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute body/local requires atom style body
Self-explanatory.
E: Invalid index in compute body/local command
Self-explanatory.
E: Invalid index for non-body particles in compute body/local command
Only indices 1,2,3 can be used for non-body particles.
*/

40
src/BODY/compute_temp_body.h
View File

@ -50,43 +50,3 @@ class ComputeTempBody : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute temp/body requires atom style body
Self-explanatory.
E: Compute temp/body requires bodies
This compute can only be applied to body particles.
E: Could not find compute ID for temperature bias
Self-explanatory.
E: Bias compute does not calculate temperature
The specified compute must compute temperature.
E: Bias compute does not calculate a velocity bias
The specified compute must compute a bias for temperature.
E: Bias compute group does not match compute group
The specified compute must operate on the same group as the parent
compute.
E: Temperature compute degrees of freedom < 0
This should not happen if you are calculating the temperature
on a valid set of atoms.
*/

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