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

Merge remote-tracking branch 'github/develop' into atom-style-var-with-python

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-05-01 16:03:27 -04:00
parent 9b36c58eb2 165f34bc12
commit e474bfeece
116 changed files with 6650 additions and 3336 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
README
View File

@ -34,6 +34,7 @@ lib additional provided or external libraries
potentials interatomic potential files
python Python module for LAMMPS
src source files
third_party Copies of thirdparty software bundled with LAMMPS
tools pre- and post-processing tools
unittest test programs for use with CTest
.github Git and GitHub related files and tools

4
cmake/CMakeLists.txt
View File

@ -44,6 +44,7 @@ set(LAMMPS_DOC_DIR ${LAMMPS_DIR}/doc)
set(LAMMPS_TOOLS_DIR ${LAMMPS_DIR}/tools)
set(LAMMPS_PYTHON_DIR ${LAMMPS_DIR}/python)
set(LAMMPS_POTENTIALS_DIR ${LAMMPS_DIR}/potentials)
set(LAMMPS_THIRDPARTY_DIR ${LAMMPS_DIR}/third_party)
set(LAMMPS_DOWNLOADS_URL "https://download.lammps.org" CACHE STRING "Base URL for LAMMPS downloads")
set(LAMMPS_POTENTIALS_URL "${LAMMPS_DOWNLOADS_URL}/potentials")
@ -131,7 +132,7 @@ endif()
# silence nvcc warnings
if((PKG_KOKKOS) AND (Kokkos_ENABLE_CUDA) AND NOT (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}" "-Xcudafe --diag_suppress=unrecognized_pragma,--diag_suppress=128")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xcudafe --diag_suppress=unrecognized_pragma,--diag_suppress=128")
endif()
# we *require* C++11 without extensions but prefer C++17.
@ -370,6 +371,7 @@ endforeach()
# packages with special compiler needs or external libs
######################################################
target_include_directories(lammps PUBLIC $<BUILD_INTERFACE:${LAMMPS_SOURCE_DIR}>)
target_include_directories(lammps PUBLIC $<BUILD_INTERFACE:${LAMMPS_THIRDPARTY_DIR}>)
if(PKG_ADIOS)
# The search for ADIOS2 must come before MPI because

2
cmake/Modules/LAMMPSUtils.cmake
View File

@ -30,7 +30,7 @@ function(check_omp_h_include)
if(OpenMP_CXX_FOUND)
set(CMAKE_REQUIRED_FLAGS ${OpenMP_CXX_FLAGS})
set(CMAKE_REQUIRED_INCLUDES ${OpenMP_CXX_INCLUDE_DIRS})
set(CMAKE_REQUIRED_LINK_OPTIONS ${OpenMP_CXX_FLAGS})
separate_arguments(CMAKE_REQUIRED_LINK_OPTIONS NATIVE_COMMAND ${OpenMP_CXX_FLAGS}) # needs to be a list
set(CMAKE_REQUIRED_LIBRARIES ${OpenMP_CXX_LIBRARIES})
# there are all kinds of problems with finding omp.h
# for Clang and derived compilers so we pretend it is there.

15
cmake/Modules/Packages/SCAFACOS.cmake
View File

@ -14,27 +14,16 @@ endif()
option(DOWNLOAD_SCAFACOS "Download ScaFaCoS library instead of using an already installed one" ${DOWNLOAD_SCAFACOS_DEFAULT})
if(DOWNLOAD_SCAFACOS)
message(STATUS "ScaFaCoS download requested - we will build our own")
set(SCAFACOS_URL "https://github.com/scafacos/scafacos/releases/download/v1.0.1/scafacos-1.0.1.tar.gz" CACHE STRING "URL for SCAFACOS tarball")
set(SCAFACOS_MD5 "bd46d74e3296bd8a444d731bb10c1738" CACHE STRING "MD5 checksum of SCAFACOS tarball")
set(SCAFACOS_URL "https://github.com/scafacos/scafacos/releases/download/v1.0.4/scafacos-1.0.4.tar.gz" CACHE STRING "URL for SCAFACOS tarball")
set(SCAFACOS_MD5 "23867540ec32e63ce71d6ecc105278d2" CACHE STRING "MD5 checksum of SCAFACOS tarball")
mark_as_advanced(SCAFACOS_URL)
mark_as_advanced(SCAFACOS_MD5)
GetFallbackURL(SCAFACOS_URL SCAFACOS_FALLBACK)
# version 1.0.1 needs a patch to compile and linke cleanly with GCC 10 and later.
file(DOWNLOAD ${LAMMPS_THIRDPARTY_URL}/scafacos-1.0.1-fix.diff ${CMAKE_CURRENT_BINARY_DIR}/scafacos-1.0.1.fix.diff
EXPECTED_HASH MD5=4baa1333bb28fcce102d505e1992d032)
find_program(HAVE_PATCH patch)
if(NOT HAVE_PATCH)
message(FATAL_ERROR "The 'patch' program is required to build the ScaFaCoS library")
endif()
include(ExternalProject)
ExternalProject_Add(scafacos_build
URL ${SCAFACOS_URL} ${SCAFACOS_FALLBACK}
URL_MD5 ${SCAFACOS_MD5}
PATCH_COMMAND patch -p1 < ${CMAKE_CURRENT_BINARY_DIR}/scafacos-1.0.1.fix.diff
CONFIGURE_COMMAND <SOURCE_DIR>/configure --prefix=<INSTALL_DIR> --disable-doc
--enable-fcs-solvers=fmm,p2nfft,direct,ewald,p3m
--with-internal-fftw --with-internal-pfft

2
cmake/presets/kokkos-sycl-intel.cmake
View File

@ -23,7 +23,7 @@ set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
# set(_intel_sycl_flags " -w -fsycl -flink-huge-device-code -fsycl-targets=spir64_gen "
set(_intel_sycl_flags " -w -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=spir64_gen "
set(_intel_sycl_flags " -w -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=spir64_gen ")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${_intel_sycl_flags}" CACHE STRING "" FORCE)
#set(CMAKE_EXE_LINKER_FLAGS "-fsycl -flink-huge-device-code -fsycl-targets=spir64_gen " CACHE STRING "" FORCE)

6
doc/graphviz/lammps-releases.dot
View File

@ -5,13 +5,13 @@ digraph releases {
github -> develop [label="Merge commits"];
{
rank = "same";
work [shape="none" label="Development branches:"]
work [shape="none" label="Development branches:" fontname="bold"]
develop [label="'develop' branch" height=0.75];
maintenance [label="'maintenance' branch" height=0.75];
};
{
rank = "same";
upload [shape="none" label="Release branches:"]
upload [shape="none" label="Release branches:" fontname="bold"]
release [label="'release' branch" height=0.75];
stable [label="'stable' branch" height=0.75];
};
@ -22,7 +22,7 @@ digraph releases {
maintenance -> stable [label="Updates to stable release"];
{
rank = "same";
tag [shape="none" label="Applied tags:"];
tag [shape="none" label="Applied tags:" fontname="bold"];
patchtag [shape="box" label="patch_<date>"];
stabletag [shape="box" label="stable_<date>"];
updatetag [shape="box" label="stable_<date>_update<num>"];

24
doc/src/Build.rst
View File

@ -14,32 +14,10 @@ As an alternative, you can download a package with pre-built executables
or automated build trees, as described in the :doc:`Install <Install>`
section of the manual.
Prerequisites
-------------
Which software you need to compile and use LAMMPS strongly depends on
which :doc:`features and settings <Build_settings>` and which
:doc:`optional packages <Packages_list>` you are trying to include.
Common to all is that you need a C++ and C compiler, where the C++
compiler has to support at least the C++11 standard (note that some
compilers require command-line flag to activate C++11 support).
Furthermore, if you are building with CMake, you need at least CMake
version 3.20 and a compatible build tool (make or ninja-build); if you
are building the the legacy GNU make based build system you need GNU
make (other make variants are not going to work since the build system
uses features unique to GNU make) and a Unix-like build environment with
a Bourne shell, and shell tools like "sed", "grep", "touch", "test",
"tr", "cp", "mv", "rm", "ln", "diff" and so on. Parts of LAMMPS
interface with or use Python version 3.6 or later.
The LAMMPS developers aim to keep LAMMPS very portable and usable -
at least in parts - on most operating systems commonly used for
running MD simulations. Please see the :doc:`section on portablility
<Intro_portability>` for more details.
.. toctree::
:maxdepth: 1
Build_prerequisites
Build_cmake
Build_make
Build_link

22
doc/src/Build_prerequisites.rst Normal file
View File

@ -0,0 +1,22 @@
Prerequisites
-------------
Which software you need to compile and use LAMMPS strongly depends on
which :doc:`features and settings <Build_settings>` and which
:doc:`optional packages <Packages_list>` you are trying to include.
Common to all is that you need a C++ and C compiler, where the C++
compiler has to support at least the C++11 standard (note that some
compilers require command-line flag to activate C++11 support).
Furthermore, if you are building with CMake, you need at least CMake
version 3.20 and a compatible build tool (make or ninja-build); if you
are building the the legacy GNU make based build system you need GNU
make (other make variants are not going to work since the build system
uses features unique to GNU make) and a Unix-like build environment with
a Bourne shell, and shell tools like "sed", "grep", "touch", "test",
"tr", "cp", "mv", "rm", "ln", "diff" and so on. Parts of LAMMPS
interface with or use Python version 3.6 or later.
The LAMMPS developers aim to keep LAMMPS very portable and usable -
at least in parts - on most operating systems commonly used for
running MD simulations. Please see the :doc:`section on portablility
<Intro_portability>` for more details.

1
doc/src/Commands_pair.rst
View File

@ -179,6 +179,7 @@ OPT.
* :doc:`lj/long/dipole/long <pair_dipole>`
* :doc:`lj/long/tip4p/long (o) <pair_lj_long>`
* :doc:`lj/mdf <pair_mdf>`
* :doc:`lj/pirani (o) <pair_lj_pirani>`
* :doc:`lj/relres (o) <pair_lj_relres>`
* :doc:`lj/spica (gko) <pair_spica>`
* :doc:`lj/spica/coul/long (gko) <pair_spica>`

6
doc/src/Commands_removed.rst
View File

@ -1,7 +1,7 @@
Removed commands and packages
=============================
.. contents:: \
.. contents::
------
@ -15,7 +15,7 @@ with the direct alternative (if available) and print a warning.
LAMMPS shell
------------
.. versionchanged:: 29Aug2024
.. deprecated:: 29Aug2024
The LAMMPS shell has been removed from the LAMMPS distribution. Users
are encouraged to use the :ref:`LAMMPS-GUI <lammps_gui>` tool instead.
@ -23,7 +23,7 @@ are encouraged to use the :ref:`LAMMPS-GUI <lammps_gui>` tool instead.
i-PI tool
---------
.. versionchanged:: 27Jun2024
.. deprecated:: 27Jun2024
The i-PI tool has been removed from the LAMMPS distribution. Instead,
instructions to install i-PI from PyPI via pip are provided.

18
doc/src/Errors_details.rst
View File

@ -159,13 +159,17 @@ angle, dihedral, or improper with just one atom in the actual
sub-domain. Typically, this cutoff is set to the largest cutoff from
the :doc:`pair style(s) <pair_style>` plus the :doc:`neighbor list skin
distance <neighbor>` and will typically be sufficient for all bonded
interactions. But if the pair style cutoff is small, this may not be
enough. LAMMPS will print a warning in this case using some heuristic
based on the equilibrium bond length, but that still may not be
sufficient for cases where the force constants are small and thus bonds
may be stretched very far. The communication cutoff can be adjusted
with :doc:`comm_modify cutoff \<value\> <comm_modify>`, but setting this
too large will waste CPU time and memory.
interactions. But if the pair style cutoff is small (e.g. with a
repulsive-only Lennard-Jones potential) this may not be enough. It is
even worse if there is no pair style defined (or the pair style is set
to "none"), since then there will be no ghost atoms created at all.
The communication cutoff can be set or adjusted with :doc:`comm_modify
cutoff \<value\> <comm_modify>`, but setting this too large will waste
CPU time and memory. LAMMPS will print warnings in these cases. For
bonds it uses some heuristic based on the equilibrium bond length, but
that still may not be sufficient for cases where the force constants are
small and thus bonds may be stretched very far.
.. _hint09:

5
doc/src/Howto_lammps_gui.rst
View File

@ -308,7 +308,10 @@ of the *Output* window showing how many warnings and errors were
detected and how many lines the entire output has. By clicking on the
button on the right with the warning symbol or by using the keyboard
shortcut `Ctrl-N` (`Command-N` on macOS), you can jump to the next
line with a warning or error.
line with a warning or error. If there is a URL pointing to additional
explanations in the online manual, that URL will be highlighted and
double-clicking on it shall open the corresponding manual page in
the web browser. The option is also available from the context menu.
By default, the *Output* window is replaced each time a run is started.
The runs are counted and the run number for the current run is displayed

13
doc/src/Install_git.rst
View File

@ -12,19 +12,10 @@ several advantages:
LAMMPS. For that, you should first create your own :doc:`fork on
GitHub <Howto_github>`, though.
You must have `git <git_>`_ installed on your system to use the
commands explained below to communicate with the git servers on
GitHub. For people still using subversion (svn), GitHub also
provides `limited support for subversion clients <svn_>`_.
.. note::
As of October 2016, the official home of public LAMMPS development is
on GitHub. The previously advertised LAMMPS git repositories on
git.lammps.org and bitbucket.org are now offline or deprecated.
You must have `git <git_>`_ installed on your system to use the commands
explained below to communicate with the git servers on GitHub.
.. _git: https://git-scm.com
.. _svn: https://help.github.com/en/github/importing-your-projects-to-github/working-with-subversion-on-github
You can follow the LAMMPS development on 4 different git branches:

BIN
doc/src/JPG/lammps-releases.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 68 KiB

After

Width:  |  Height:  |  Size: 70 KiB

4
doc/src/Modify_pair.rst
View File

@ -46,6 +46,8 @@ Here is a brief list of some the class methods in the Pair class that
+---------------------------------+------------------------------------------------------------------------+
| compute_inner/middle/outer | versions of compute used by rRESPA |
+---------------------------------+------------------------------------------------------------------------+
| compute_atomic_energy | energy of one atom, equivalent to per-atom energy |
+---------------------------------+------------------------------------------------------------------------+
| memory_usage | return estimated amount of memory used by the pair style |
+---------------------------------+------------------------------------------------------------------------+
| modify_params | process arguments to pair_modify command |
@ -122,3 +124,5 @@ setting.
+---------------------------------+-------------------------------------------------------------+---------+
| spinflag | 1 if compatible with spin kspace_style | 0 |
+---------------------------------+-------------------------------------------------------------+---------+
| atomic_energy_enable | 1 if compute_atomic_energy() routine exists | 0 |
+---------------------------------+-------------------------------------------------------------+---------+

10
doc/src/Tools.rst
View File

@ -1250,10 +1250,10 @@ tabulate tool
.. versionadded:: 22Dec2022
The ``tabulate`` folder contains Python scripts scripts to generate tabulated
potential files for LAMMPS. The bulk of the code is in the ``tabulate`` module
in the ``tabulate.py`` file. Some example files demonstrating its use are
included. See the README file for more information.
The ``tabulate`` folder contains Python scripts scripts to generate and
visualize tabulated potential files for LAMMPS. The bulk of the code is in the
``tabulate`` module in the ``tabulate.py`` file. Some example files
demonstrating its use are included. See the README file for more information.
----------
@ -1276,7 +1276,7 @@ Those scripts were written by Steve Plimpton sjplimp at gmail.com
valgrind tool
-------------
The ``valgrind`` folder contains additional suppressions fur LAMMPS when
The ``valgrind`` folder contains additional suppressions for LAMMPS when
using `valgrind's <https://valgrind.org/>`_ ` `memcheck tool
<https://valgrind.org/info/tools.html#memcheck>`_ to search for memory
access violation and memory leaks. These suppressions are automatically

3
doc/src/bond_bpm_rotational.rst
View File

@ -215,6 +215,9 @@ for an overview of LAMMPS output options.
The vector or array will be floating point values that correspond to
the specified attribute.
Any settings with the *store/local* option are not saved to a restart
file and must be redefined.
The single() function of this bond style returns 0.0 for the energy
of a bonded interaction, since energy is not conserved in these
dissipative potentials. It also returns only the normal component of

3
doc/src/bond_bpm_spring.rst
View File

@ -215,6 +215,9 @@ for an overview of LAMMPS output options.
The vector or array will be floating point values that correspond to
the specified attribute.
Any settings with the *store/local* option are not saved to a restart
file and must be redefined.
The potential energy and the single() function of this bond style return
:math:`k (r - r_0)^2 / 2` as a proxy of the energy of a bonded interaction,
ignoring any volumetric/smoothing factors or dissipative forces. The single()

2
doc/src/compute_rheo_property_atom.rst
View File

@ -88,7 +88,7 @@ The *phase* property indicates whether the particle is in a fluid state,
a value of 0, or a solid state, a value of 1.
The *surface* property indicates the surface designation produced by
the *interface/reconstruct* option of :doc:`fix rheo <fix_rheo>`. Bulk
the *surface/detection* option of :doc:`fix rheo <fix_rheo>`. Bulk
particles have a value of 0, surface particles have a value of 1, and
splash particles have a value of 2. The *surface/r* property is the
distance from the surface, up to the kernel cutoff length. Surface particles

2
doc/src/fix_adapt.rst
View File

@ -214,6 +214,8 @@ formulas for the meaning of these parameters:
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`lj/mdf <pair_mdf>` | epsilon,sigma | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`lj/pirani <pair_lj_pirani>` | alpha, beta, gamma, rm, epsilon | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`lj/sf/dipole/sf <pair_dipole>` | epsilon,sigma,scale | type pairs |
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
| :doc:`lubricate <pair_lubricate>` | mu | global |

43
doc/src/fix_sgcmc.rst
View File

@ -30,7 +30,9 @@ Syntax
N = number of times sampling window is moved during one MC cycle
*window_size* frac
frac = size of sampling window (must be between 0.5 and 1.0)
*atomic/energy* yes/no
yes = use the atomic energy method to calculate energy changes
no = use the default method to calculate energy changes
Examples
""""""""
@ -127,6 +129,14 @@ The number of times the window is moved during a MC cycle is set using
the parameter *window_moves* (see Sect. III.B in :ref:`Sadigh1
<Sadigh1>` for details).
The *atomic/energy* keyword controls which method is used for calculating
the energy change when atom types are swapped. A value of *no*
uses the default method, see discussion below in Restrictions section.
A value of *yes* uses the atomic energy method,
if the method has been implemented for the LAMMPS energy model,
otherwise LAMMPS will exit with an error message.
So far this has only been implemented for EAM type potentials.
------------
Restart, fix_modify, output, run start/stop, minimize info
@ -159,16 +169,26 @@ page for more info.
This fix style requires an :doc:`atom style <atom_style>` with per atom
type masses.
At present the fix provides optimized subroutines for EAM type
potentials (see above) that calculate potential energy changes due to
*local* atom type swaps very efficiently. Other potentials are
supported by using the generic potential functions. This, however, will
lead to exceedingly slow simulations since it implies that the
energy of the *entire* system is recomputed at each MC trial step. If
other potentials are to be used it is strongly recommended to modify and
optimize the existing generic potential functions for this purpose.
Also, the generic energy calculation can not be used for parallel
execution i.e. it only works with a single MPI process.
The fix provides three methods for calculating the potential energy
change due to atom type swaps. For EAM type potentials, the default
method is a carefully optimized local energy change calculation that
is part of the source code for this fix. It takes advantage of the
specific computational and communication requirements of EAM. Customizing
the local method to handle other energy models such as Tersoff has been done,
but these cases are not supported in the public LAMMPS code.
For all other LAMMPS energy models, the default method calculates
the *total* potential energy of the system before and after each
atom type swap. This method does not depend on the details of the
energy model and so is guaranteed to be correct. It is also
orders of magnitude slower than the custom EAM calculation.
In addition, it can not be used with parallel execution i.e. only
a single MPI process is allowed.
The third method uses the *atomic/energy* keyword described above.
This allows parallel execution and it is also a local calculation,
making it only a bit slower than a fully-optimized local calculation.
So far, this has been implemented for EAM type potentials.
It is straightforward to extend this to other potentials,
requiring adding an atomic energy method to the pair style.
------------
@ -180,6 +200,7 @@ The optional parameters default to the following values:
* *randseed* = 324234
* *window_moves* = 8
* *window_size* = automatic
* *atomic/energy* = no
------------

163
doc/src/pair_lj_pirani.rst Normal file
View File

@ -0,0 +1,163 @@
.. index:: pair_style lj/pirani
.. index:: pair_style lj/pirani/omp
pair_style lj/pirani command
============================
Accelerator Variants: *lj/pirani/omp*
Syntax
""""""
.. code-block:: LAMMPS
pair_style lj/pirani cutoff
* lj/pirani = name of the pair style
* cutoff = global cutoff (distance units)
Examples
""""""""
.. code-block:: LAMMPS
pair_style lj/pirani 10.0
pair_coeff 1 1 4.0 7.0 6.0 3.5 0.0045
Description
"""""""""""
.. versionadded:: TBD
Pair style *lj/pirani* computes pairwise interactions from an Improved
Lennard-Jones (ILJ) potential according to :ref:`(Pirani) <Pirani>`.
The ILJ force field is adequate to model both equilibrium and
non-equilibrium properties of matter, in gaseous and condensed phases,
and at gas-surface interfaces. In particular, its use improves the
description of elementary process dynamics where the traditional
Lennard-Jones (LJ) formulation is usually applied.
.. math::
x = r/R_m \\
n_x = \alpha*x^2 + \beta \\
\gamma \equiv m \\
V(x) = \varepsilon \cdot \left( \frac{\gamma}{ n_x - \gamma} \left(\frac{1}{x} \right)^{n_x}
- \frac{n_x}{n_x - \gamma} \left(\frac{1}{x} \right)^{\gamma} \right) \qquad r < r_c
:math:`r_c` is the cutoff.
An additional parameter, :math:`\alpha`, has been introduced in order to
be able to recover the traditional Lennard-Jones 12-6 with a specific
choice of parameters. With :math:`R_m \equiv r_0 = \sigma \cdot 2^{1 /
6}`, :math:`\alpha = 0`, :math:`\beta = 12` and :math:`\gamma = 6` it is
straightforward to prove that LJ 12-6 is obtained. Also, it can be
verified that using :math:`\alpha= 4`, :math:`\beta= 8` and
:math:`\gamma = 6`, at the equilibrium distance, the first and second
derivatives of ILJ match those of LJ 12-6. The parameter :math:`R_m`
corresponds to the equilibrium distance and :math:`\epsilon` to the well
depth.
This potential provides some advantages with respect to the standard LJ
potential, as explained in :ref:`(Pirani) <Pirani>`: it provides a more
realistic description of the long range behavior and an attenuation of
the hardness of the repulsive wall.
This force field can be used for neutral-neutral (:math:`\gamma = 6`),
ion-neutral (:math:`\gamma = 4`) or ion-ion systems (:math:`\gamma =
1`). Notice that this implementation does not include explicit
electrostatic interactions. If these are desired, this pair style
should be used along with a Coulomb pair style like
:doc:`pair styles coul/cut or coul/long <pair_coul>` by using
:doc:`pair style hybrid/overlay <pair_hybrid>` and a suitable
:doc:`kspace style <kspace_style>`, if needed.
As discussed in :ref:`(Pirani) <Pirani>`, analysis of a variety of
systems showed that :math:`\alpha= 4` generally works very well. In
some special cases (e.g. those involving very small multiple charged
ions) this factor may take a slightly different value. The parameter
:math:`\beta` codifies the hardness (polarizability) of the interacting
partners, and for neutral-neutral systems it usually ranges from 6
to 11. Moreover, the modulation of :math:`\beta` can model additional
interaction effects, such as charge transfer in the perturbative limit,
and can mitigate the effect of some uncertainty in the data used to
build up the potential function.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
above, or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands:
* :math:`\alpha` (dimensionless)
* :math:`\beta` (dimensionless)
* :math:`\gamma` (dimensionless)
* :math:`R_m` (distance units)
* :math:`\epsilon` (energy units)
* cutoff (distance units)
The last coefficient is optional. If not specified, the global cutoff is used.
----------
.. include:: accel_styles.rst
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair style does not support mixing. Thus, coefficients for all I,J
pairs must be specified explicitly.
This pair style supports the :doc:`pair_modify <pair_modify>` shift
option for the energy of the pair interaction.
The :doc:`pair_modify <pair_modify>` table options are not relevant for
this pair style.
This pair style does not support the :doc:`pair_modify <pair_modify>`
tail option for adding long-range tail corrections to energy and
pressure.
This pair style writes its information to :doc:`binary restart files
<restart>`, so pair_style and pair_coeff commands do not need to be
specified in an input script that reads a restart file.
This pair style supports the use of the *inner*, *middle*, and
*outer* keywords of the :doc:`run_style respa <run_style>` command,
meaning the pairwise forces can be partitioned by distance at different
levels of the rRESPA hierarchy. See the :doc:`run_style <run_style>`
command for details.
----------
Restrictions
""""""""""""
This pair style is only enabled if LAMMPS was built with the EXTRA-PAIR
package. See the :doc:`Build package <Build_package>` page for more
info.
Related commands
""""""""""""""""
* :doc:`pair_coeff <pair_coeff>`
* :doc:`pair_style lj/cut <pair_lj>`
Default
"""""""
none
--------------
.. _Pirani:
**(Pirani)** F. Pirani, S. Brizi, L. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi,
Phys. Chem. Chem. Phys., 2008, 10, 5489-5503.

1
doc/src/pair_style.rst
View File

@ -272,6 +272,7 @@ accelerated styles exist.
* :doc:`lj/long/dipole/long <pair_dipole>` - long-range LJ and long-range point dipoles
* :doc:`lj/long/tip4p/long <pair_lj_long>` - long-range LJ and long-range Coulomb for TIP4P water
* :doc:`lj/mdf <pair_mdf>` - LJ potential with a taper function
* :doc:`lj/pirani <pair_lj_pirani>` - Improved LJ potential
* :doc:`lj/relres <pair_lj_relres>` - LJ using multiscale Relative Resolution (RelRes) methodology :ref:`(Chaimovich) <Chaimovich2>`.
* :doc:`lj/spica <pair_spica>` - LJ for SPICA coarse-graining
* :doc:`lj/spica/coul/long <pair_spica>` - LJ for SPICA coarse-graining with long-range Coulomb

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

@ -389,6 +389,7 @@ Bretonnet
Briels
Brien
Brilliantov
Brizi
Broadwell
Broglie
brownian
@ -431,10 +432,12 @@ Camiloni
Campana
Cangi
Cao
Cappelletti
Capolungo
Caro
cartesian
Cas
Casavecchia
CasP
Caswell
Cates
@ -2922,6 +2925,7 @@ perp
Perram
persp
Persp
perturbative
peru
Peskin
Pettifor
@ -2960,6 +2964,8 @@ pimdb
Piola
pIp
pipelining
Pirani
pirani
Pisarev
Pishevar
Pitera
@ -3347,6 +3353,7 @@ Rockett
rocksalt
Rodrigues
Rohart
Roncaratti
Ronchetti
Ronevich
Rosati
@ -4060,6 +4067,7 @@ vdW
vdwl
vec
Vecchio
Vecchiocattivi
vectorial
vectorization
Vectorization

1
examples/COUPLE/plugin/liblammpsplugin.c
View File

@ -144,6 +144,7 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib)
ADDSYM(find_pair_neighlist);
ADDSYM(find_fix_neighlist);
ADDSYM(find_compute_neighlist);
ADDSYM(request_single_neighlist);
ADDSYM(neighlist_num_elements);
ADDSYM(neighlist_element_neighbors);

18
examples/COUPLE/plugin/liblammpsplugin.h
View File

@ -94,6 +94,17 @@ enum _LMP_VAR_CONST {
LMP_VAR_STRING = 3 /*!< return value will be a string (catch-all) */
};
/** Neighbor list settings constants
*
* Must be kept in sync with the equivalent constants in ``python/lammps/constants.py``,
* ``fortran/lammps.f90``, ``tools/swig/lammps.i``, and
* ``examples/COUPLE/plugin/liblammpsplugin.h`` */
enum _LMP_NEIGH_CONST {
LMP_NEIGH_HALF = 0, /*!< request (default) half neighbor list */
LMP_NEIGH_FULL = 1, /*!< request full neighbor list */
};
#ifdef __cplusplus
extern "C" {
#endif
@ -189,14 +200,17 @@ struct _liblammpsplugin {
* caller must match to how LAMMPS library is built */
#if !defined(LAMMPS_BIGBIG)
int (*create_atoms)(void *, int, int *, int *, double *, double *, int *, int);
int (*create_atoms)(void *, int, const int *, const int *, const double *, const double *,
const int *, int);
#else
int (*create_atoms)(void *, int, int64_t *, int *, double *, double *, int64_t *, int);
int (*create_atoms)(void *, int, const int64_t *, const int *, const double *, const double *,
const int64_t *, int);
#endif
int (*find_pair_neighlist)(void *, const char *, int, int, int);
int (*find_fix_neighlist)(void *, const char *, int);
int (*find_compute_neighlist)(void *, const char *, int);
int (*request_single_neighlist)(void *, const char *, int, double);
int (*neighlist_num_elements)(void *, int);
void (*neighlist_element_neighbors)(void *, int, int, int *, int *, int **);

18
examples/granular/data.particles
View File

@ -1,18 +0,0 @@
Python generated LAMMPS data file
2 atoms
1 atom types
0 0.08 xlo xhi
0 0.04 ylo yhi
0 0.08 zlo zhi
Atoms # sphere
1 1 0.004 2500 0.04 0.02 0.04 0 0 0
2 1 0.004 2500 0.04 0.02 0.04416 0 0 0
Velocities
1 0.0 0.0 1 0 0 0
2 0.0 0.0 -1 0 0 0

47
examples/granular/in.restitution
View File

@ -1,24 +1,43 @@
units si
atom_style sphere
comm_modify vel yes
boundary p p p
boundary p p f
region box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4
region box block -0.01 0.01 -0.01 0.01 0 0.08
create_box 2 box
read_data data.particles add append
group mb type 1
create_atoms 1 single 0.0 0.0 0.04
create_atoms 1 single 0.0 0.0 0.04416
set group all diameter 0.004 density 2500
pair_style granular
pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution
# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping coeff_restitution
comm_modify vel yes
pair_style granular
pair_coeff * * hertz/material 1e6 0.8 0.4 &
tangential mindlin NULL 0.0 0.5 &
damping coeff_restitution
#pair_coeff * * hooke 1e6 0.5 &
# tangential mindlin 1 1.0 0.0 &
# damping coeff_restitution
timestep 1e-9
fix 1 all nve/sphere
compute s all stress/atom NULL pair
#dump 1 all custom 2000000 op.dump id x y z vx vy vz
#dump_modify 1 pad 8
thermo_style custom step ke
run_style verlet
run 10000000
group a1 id 1
group a2 id 2
velocity a1 set 0 0 1.0
velocity a2 set 0 0 -1.0
compute z1 a1 reduce sum z
compute z2 a2 reduce sum z
compute v1 a1 reduce sum vz
compute v2 a2 reduce sum vz
compute f1 a1 reduce sum fz
compute f2 a2 reduce sum fz
variable dz equal c_z1 - c_z2
# dump 1 all custom 2000000 op.dump id x y z vx vy vz
thermo 10000
thermo_style custom step ke v_dz c_v1 c_v2 c_f1 c_f2
run 1000000

80
examples/granular/log.13May23.restitution.g++.1
View File

@ -1,80 +0,0 @@
LAMMPS (17 Apr 2024 - Development - patch_17Apr2024-93-g4e7bddaa0b)
using 1 OpenMP thread(s) per MPI task
units si
atom_style sphere
boundary p p f
region box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4
create_box 2 box
Created orthogonal box = (0 0 0) to (0.08 0.04 0.08)
1 by 1 by 1 MPI processor grid
read_data data.particles add append
Reading data file ...
orthogonal box = (0 0 0) to (0.08 0.04 0.08)
1 by 1 by 1 MPI processor grid
reading atoms ...
2 atoms
reading velocities ...
2 velocities
read_data CPU = 0.002 seconds
group mb type 1
2 atoms in group mb
pair_style granular
pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution
# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping coeff_restitution
comm_modify vel yes
timestep 1e-9
fix 1 all nve/sphere
compute s all stress/atom NULL pair
#dump 1 all custom 2000000 op.dump id x y z vx vy vz
#dump_modify 1 pad 8
thermo_style custom step ke
run_style verlet
run 10000000
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 0.005
ghost atom cutoff = 0.005
binsize = 0.0025, bins = 32 16 32
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair granular, perpetual
attributes: half, newton on, size, history
pair build: half/size/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 10.1 | 10.1 | 10.1 Mbytes
Step KinEng
0 8.3775804e-05
10000000 5.3616513e-05
Loop time of 5.99782 on 1 procs for 10000000 steps with 2 atoms
77.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.60235 | 0.60235 | 0.60235 | 0.0 | 10.04
Neigh | 0.00021965 | 0.00021965 | 0.00021965 | 0.0 | 0.00
Comm | 1.7939 | 1.7939 | 1.7939 | 0.0 | 29.91
Output | 2.5955e-05 | 2.5955e-05 | 2.5955e-05 | 0.0 | 0.00
Modify | 1.7622 | 1.7622 | 1.7622 | 0.0 | 29.38
Other | | 1.839 | | | 30.66
Nlocal: 2 ave 2 max 2 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 0
Ave neighs/atom = 0
Neighbor list builds = 14
Dangerous builds = 0
Total wall time: 0:00:06

195
examples/granular/log.4Feb25.restitution.g++.1 Normal file
View File

@ -0,0 +1,195 @@
LAMMPS (4 Feb 2025)
units si
atom_style sphere
comm_modify vel yes
boundary p p p
region box block -0.01 0.01 -0.01 0.01 0 0.08
create_box 2 box
Created orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08)
1 by 1 by 1 MPI processor grid
create_atoms 1 single 0.0 0.0 0.04
Created 1 atoms
using lattice units in orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08)
create_atoms CPU = 0.000 seconds
create_atoms 1 single 0.0 0.0 0.04416
Created 1 atoms
using lattice units in orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08)
create_atoms CPU = 0.000 seconds
set group all diameter 0.004 density 2500
Setting atom values ...
2 settings made for diameter
2 settings made for density
pair_style granular
pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution
#pair_coeff * * hooke 1e6 0.5 # tangential mindlin 1 1.0 0.0 # damping coeff_restitution
timestep 1e-9
fix 1 all nve/sphere
group a1 id 1
1 atoms in group a1
group a2 id 2
1 atoms in group a2
velocity a1 set 0 0 1.0
velocity a2 set 0 0 -1.0
compute z1 a1 reduce sum z
compute z2 a2 reduce sum z
compute v1 a1 reduce sum vz
compute v2 a2 reduce sum vz
compute f1 a1 reduce sum fz
compute f2 a2 reduce sum fz
variable dz equal c_z1 - c_z2
# dump 1 all custom 2000000 op.dump id x y z vx vy vz
thermo 10000
thermo_style custom step ke v_dz c_v1 c_v2 c_f1 c_f2
run 1000000
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 0.005
ghost atom cutoff = 0.005
binsize = 0.0025, bins = 8 8 32
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair granular, perpetual
attributes: half, newton on, size, history
pair build: half/size/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 9.998 | 9.998 | 9.998 Mbytes
Step KinEng v_dz c_v1 c_v2 c_f1 c_f2
0 8.3775804e-05 -0.00416 1 -1 0 0
10000 8.3775804e-05 -0.00414 1 -1 0 0
20000 8.3775804e-05 -0.00412 1 -1 0 0
30000 8.3775804e-05 -0.0041 1 -1 0 0
40000 8.3775804e-05 -0.00408 1 -1 0 0
50000 8.3775804e-05 -0.00406 1 -1 0 0
60000 8.3775804e-05 -0.00404 1 -1 0 0
70000 8.3775804e-05 -0.00402 1 -1 0 0
80000 8.3775804e-05 -0.004 1 -1 0 0
90000 8.3411065e-05 -0.003980019 0.99782075 -0.99782075 -0.023914747 0.023914747
100000 8.2852688e-05 -0.0039600945 0.9944753 -0.9944753 -0.032005131 0.032005131
110000 8.2139641e-05 -0.0039402463 0.99018672 -0.99018672 -0.039875404 0.039875404
120000 8.1272296e-05 -0.0039204934 0.98494496 -0.98494496 -0.048007824 0.048007824
130000 8.0246788e-05 -0.0039008551 0.97871113 -0.97871113 -0.056503872 0.056503872
140000 7.9058986e-05 -0.0038813518 0.97144075 -0.97144075 -0.065373554 0.065373554
150000 7.7705654e-05 -0.0038620047 0.9630903 -0.9630903 -0.074595879 0.074595879
160000 7.6184906e-05 -0.0038428357 0.95361959 -0.95361959 -0.084137355 0.084137355
170000 7.4496418e-05 -0.0038238676 0.94299284 -0.94299284 -0.093958893 0.093958893
180000 7.2641536e-05 -0.0038051239 0.93117907 -0.93117907 -0.10401872 0.10401872
190000 7.0623328e-05 -0.0037866285 0.91815243 -0.91815243 -0.11427372 0.11427372
200000 6.8446602e-05 -0.003768406 0.90389221 -0.90389221 -0.12468011 0.12468011
210000 6.6117901e-05 -0.0037504812 0.88838298 -0.88838298 -0.13519381 0.13519381
220000 6.3645478e-05 -0.0037328791 0.87161455 -0.87161455 -0.14577066 0.14577066
230000 6.1039243e-05 -0.003715625 0.85358204 -0.85358204 -0.1563666 0.1563666
240000 5.8310702e-05 -0.0036987442 0.83428576 -0.83428576 -0.16693776 0.16693776
250000 5.5472871e-05 -0.003682262 0.8137313 -0.8137313 -0.17744066 0.17744066
260000 5.2540172e-05 -0.0036662033 0.79192936 -0.79192936 -0.18783225 0.18783225
270000 4.9528314e-05 -0.003650593 0.76889577 -0.76889577 -0.19807012 0.19807012
280000 4.6454158e-05 -0.0036354555 0.74465137 -0.74465137 -0.20811258 0.20811258
290000 4.3335566e-05 -0.0036208149 0.71922195 -0.71922195 -0.21791884 0.21791884
300000 4.0191232e-05 -0.0036066944 0.69263807 -0.69263807 -0.22744912 0.22744912
310000 3.7040511e-05 -0.0035931168 0.66493499 -0.66493499 -0.23666485 0.23666485
320000 3.390323e-05 -0.0035801042 0.63615249 -0.63615249 -0.24552878 0.24552878
330000 3.0799488e-05 -0.0035676776 0.60633473 -0.60633473 -0.25400513 0.25400513
340000 2.7749462e-05 -0.0035558573 0.57553002 -0.57553002 -0.26205979 0.26205979
350000 2.4773197e-05 -0.0035446626 0.54379063 -0.54379063 -0.2696604 0.2696604
360000 2.1890403e-05 -0.0035341116 0.51117262 -0.51117262 -0.27677654 0.27677654
370000 1.9120254e-05 -0.0035242212 0.47773551 -0.47773551 -0.28337983 0.28337983
380000 1.6481181e-05 -0.0035150072 0.44354212 -0.44354212 -0.28944409 0.28944409
390000 1.3990689e-05 -0.0035064841 0.40865822 -0.40865822 -0.29494544 0.29494544
400000 1.1665166e-05 -0.003498665 0.37315233 -0.37315233 -0.2998624 0.2998624
410000 9.5197195e-06 -0.0034915617 0.33709536 -0.33709536 -0.304176 0.304176
420000 7.5680136e-06 -0.0034851844 0.30056032 -0.30056032 -0.30786987 0.30786987
430000 5.8221324e-06 -0.003479542 0.26362205 -0.26362205 -0.31093026 0.31093026
440000 4.2924559e-06 -0.0034746417 0.22635684 -0.22635684 -0.31334618 0.31334618
450000 2.9875585e-06 -0.0034704894 0.18884214 -0.18884214 -0.31510936 0.31510936
460000 1.9141264e-06 -0.0034670891 0.15115621 -0.15115621 -0.31621432 0.31621432
470000 1.0768988e-06 -0.0034644437 0.11337783 -0.11337783 -0.31665837 0.31665837
480000 4.786302e-07 -0.0034625541 0.075585893 -0.075585893 -0.31644161 0.31644161
490000 1.2007709e-07 -0.0034614198 0.037859142 -0.037859142 -0.31556689 0.31556689
500000 6.3727744e-12 -0.0034610388 0.0002758068 -0.0002758068 -0.31403982 0.31403982
510000 1.1522726e-07 -0.0034614073 -0.03708671 0.03708671 -0.31186864 0.31186864
520000 4.6064472e-07 -0.0034625203 -0.074152149 0.074152149 -0.30906426 0.30906426
530000 1.029334e-06 -0.0034643709 -0.1108457 0.1108457 -0.30564009 0.30564009
540000 1.812635e-06 -0.0034669512 -0.14709431 0.14709431 -0.30161199 0.30161199
550000 2.8002645e-06 -0.0034702513 -0.18282695 0.18282695 -0.29699817 0.29699817
560000 3.9804448e-06 -0.0034742603 -0.21797491 0.21797491 -0.29181905 0.29181905
570000 5.3400475e-06 -0.0034789659 -0.25247202 0.25247202 -0.28609717 0.28609717
580000 6.8647484e-06 -0.0034843545 -0.28625495 0.28625495 -0.27985701 0.27985701
590000 8.5391931e-06 -0.003490411 -0.31926339 0.31926339 -0.2731249 0.2731249
600000 1.0347171e-05 -0.0034971194 -0.35144026 0.35144026 -0.2659288 0.2659288
610000 1.2271792e-05 -0.0035044627 -0.38273192 0.38273192 -0.25829822 0.25829822
620000 1.429567e-05 -0.0035124225 -0.41308835 0.41308835 -0.25026402 0.25026402
630000 1.6401103e-05 -0.0035209797 -0.44246327 0.44246327 -0.24185823 0.24185823
640000 1.8570255e-05 -0.0035301142 -0.47081428 0.47081428 -0.23311393 0.23311393
650000 2.078533e-05 -0.0035398052 -0.498103 0.498103 -0.22406507 0.22406507
660000 2.3028745e-05 -0.003550031 -0.52429514 0.52429514 -0.21474628 0.21474628
670000 2.5283288e-05 -0.0035607695 -0.54936056 0.54936056 -0.20519274 0.20519274
680000 2.7532278e-05 -0.0035719977 -0.57327337 0.57327337 -0.19544002 0.19544002
690000 2.9759697e-05 -0.0035836926 -0.59601193 0.59601193 -0.18552392 0.18552392
700000 3.1950329e-05 -0.0035958303 -0.61755887 0.61755887 -0.17548034 0.17548034
710000 3.408987e-05 -0.0036083869 -0.63790112 0.63790112 -0.16534511 0.16534511
720000 3.6165032e-05 -0.0036213382 -0.65702988 0.65702988 -0.15515392 0.15515392
730000 3.8163631e-05 -0.00363466 -0.67494058 0.67494058 -0.14494218 0.14494218
740000 4.0074659e-05 -0.0036483277 -0.69163285 0.69163285 -0.13474489 0.13474489
750000 4.1888343e-05 -0.0036623172 -0.7071105 0.7071105 -0.1245966 0.1245966
760000 4.3596185e-05 -0.0036766041 -0.7213814 0.7213814 -0.11453133 0.11453133
770000 4.5190996e-05 -0.0036911645 -0.73445747 0.73445747 -0.1045825 0.1045825
780000 4.6666906e-05 -0.0037059745 -0.74635458 0.74635458 -0.094782939 0.094782939
790000 4.8019368e-05 -0.0037210109 -0.75709246 0.75709246 -0.085164857 0.085164857
800000 4.9245153e-05 -0.0037362507 -0.76669467 0.76669467 -0.075759891 0.075759891
810000 5.0342325e-05 -0.0037516713 -0.77518852 0.77518852 -0.066599178 0.066599178
820000 5.1310215e-05 -0.003767251 -0.78260499 0.78260499 -0.057713474 0.057713474
830000 5.2149388e-05 -0.0037829686 -0.78897875 0.78897875 -0.049133335 0.049133335
840000 5.2861601e-05 -0.0037988035 -0.79434809 0.79434809 -0.040889384 0.040889384
850000 5.3449761e-05 -0.0038147361 -0.79875498 0.79875498 -0.03301269 0.03301269
860000 5.3917883e-05 -0.0038307476 -0.80224517 0.80224517 -0.025535302 0.025535302
870000 5.4271056e-05 -0.0038468201 -0.80486832 0.80486832 -0.018491033 0.018491033
880000 5.4515415e-05 -0.0038629369 -0.80667827 0.80667827 -0.011916591 0.011916591
890000 5.4658137e-05 -0.0038790822 -0.80773352 0.80773352 -0.00585333 0.00585333
900000 5.4707463e-05 -0.0038952416 -0.80809791 0.80809791 -0.00035001143 0.00035001143
910000 5.4672786e-05 -0.003911402 -0.80784176 0.80784176 0.0045325009 -0.0045325009
920000 5.4564826e-05 -0.0039275517 -0.80704376 0.80704376 0.0087126344 -0.0087126344
930000 5.4396015e-05 -0.0039436807 -0.80579439 0.80579439 0.012070422 -0.012070422
940000 5.4181322e-05 -0.0039597811 -0.80420264 0.80420264 0.014404618 -0.014404618
950000 5.3940325e-05 -0.0039758475 -0.80241211 0.80241211 0.015295834 -0.015295834
960000 5.3704817e-05 -0.0039918778 -0.8006585 0.8006585 0.013305187 -0.013305187
970000 5.3616513e-05 -0.0040078808 -0.79999999 0.79999999 0 0
980000 5.3616513e-05 -0.0040238808 -0.79999999 0.79999999 0 0
990000 5.3616513e-05 -0.0040398808 -0.79999999 0.79999999 0 0
1000000 5.3616513e-05 -0.0040558808 -0.79999999 0.79999999 0 0
Loop time of 0.572025 on 1 procs for 1000000 steps with 2 atoms
99.3% 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.1248 | 0.1248 | 0.1248 | 0.0 | 21.82
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.21166 | 0.21166 | 0.21166 | 0.0 | 37.00
Output | 0.00057419 | 0.00057419 | 0.00057419 | 0.0 | 0.10
Modify | 0.12695 | 0.12695 | 0.12695 | 0.0 | 22.19
Other | | 0.108 | | | 18.89
Nlocal: 2 ave 2 max 2 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1 ave 1 max 1 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 1
Ave neighs/atom = 0.5
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:01

52
fortran/lammps.f90
View File

@ -74,7 +74,9 @@ MODULE LIBLAMMPS
LMP_VAR_EQUAL = 0, & ! equal-style variables (and compatible)
LMP_VAR_ATOM = 1, & ! atom-style variables
LMP_VAR_VECTOR = 2, & ! vector variables
LMP_VAR_STRING = 3 ! string variables (everything else)
LMP_VAR_STRING = 3, & ! string variables (everything else)
LMP_NEIGH_HALF = 0, & ! request (default) half neighbor list
LMP_NEIGH_FULL = 1 ! request full neighbor list
! Constants we set once (in the constructor) and never need to check again
INTEGER(c_int), SAVE :: SIZE_TAGINT, SIZE_BIGINT, SIZE_IMAGEINT
@ -195,10 +197,11 @@ MODULE LIBLAMMPS
PROCEDURE, PRIVATE :: lmp_create_atoms_bigbig
GENERIC :: create_atoms => lmp_create_atoms_int, &
lmp_create_atoms_bigbig
PROCEDURE :: find_pair_neighlist => lmp_find_pair_neighlist
PROCEDURE :: find_fix_neighlist => lmp_find_fix_neighlist
PROCEDURE :: find_compute_neighlist => lmp_find_compute_neighlist
PROCEDURE :: neighlist_num_elements => lmp_neighlist_num_elements
PROCEDURE :: find_pair_neighlist => lmp_find_pair_neighlist
PROCEDURE :: find_fix_neighlist => lmp_find_fix_neighlist
PROCEDURE :: find_compute_neighlist => lmp_find_compute_neighlist
PROCEDURE :: request_single_neighlist => lmp_request_single_neighlist
PROCEDURE :: neighlist_num_elements => lmp_neighlist_num_elements
PROCEDURE :: neighlist_element_neighbors => lmp_neighlist_element_neighbors
PROCEDURE :: version => lmp_version
PROCEDURE, NOPASS :: get_os_info => lmp_get_os_info
@ -778,6 +781,15 @@ MODULE LIBLAMMPS
INTEGER(c_int) :: lammps_find_compute_neighlist
END FUNCTION lammps_find_compute_neighlist
FUNCTION lammps_request_single_neighlist(handle, id, flags, cutoff) BIND(C)
IMPORT :: c_int, c_double, c_ptr
IMPLICIT NONE
TYPE(c_ptr), VALUE :: handle, id
INTEGER(c_int), VALUE :: flags
REAL(c_double), VALUE :: cutoff
INTEGER(c_int) :: lammps_request_single_neighlist
END FUNCTION lammps_request_single_neighlist
FUNCTION lammps_neighlist_num_elements(handle, idx) BIND(C)
IMPORT :: c_ptr, c_int
TYPE(c_ptr), VALUE :: handle
@ -2942,6 +2954,36 @@ CONTAINS
CALL lammps_free(Cid)
END FUNCTION lmp_find_compute_neighlist
! equivalent function to lammps_request_single_neighlist
INTEGER(c_int) FUNCTION lmp_request_single_neighlist(self, id, flags, cutoff) RESULT(idx)
CLASS(lammps), INTENT(IN) :: self
CHARACTER(LEN=*), INTENT(IN) :: id
INTEGER(c_int), INTENT(IN), OPTIONAL :: flags
REAL(c_double), INTENT(IN), OPTIONAL :: cutoff
TYPE(c_ptr) :: Cid
INTEGER(c_int) :: Cflags
REAL(c_double) :: Ccutoff
IF (PRESENT(flags)) THEN
Cflags = flags
ELSE
Cflags = LMP_NEIGH_HALF
END IF
IF (PRESENT(cutoff)) THEN
Ccutoff = cutoff
ELSE
Ccutoff = 1.0_c_double
END IF
Cid = f2c_string(id)
idx = lammps_request_single_neighlist(self%handle, Cid, Cflags, Ccutoff)
IF (idx < 0) THEN
CALL lmp_error(self, LMP_ERROR_WARNING + LMP_ERROR_WORLD, &
'neighbor list build failed [Fortran/request_single_neighlist]')
END IF
CALL lammps_free(Cid)
END FUNCTION lmp_request_single_neighlist
INTEGER(c_int) FUNCTION lmp_neighlist_num_elements(self, idx) RESULT(inum)
CLASS(lammps), INTENT(IN) :: self
INTEGER(c_int), INTENT(IN) :: idx

5
lib/scafacos/Install.py
View File

@ -16,11 +16,12 @@ parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wra
# settings
version = "1.0.1"
version = "1.0.4"
# known checksums for different ScaFaCoS versions. used to validate the download.
checksums = { \
'1.0.1' : 'bd46d74e3296bd8a444d731bb10c1738' \
'1.0.1' : 'bd46d74e3296bd8a444d731bb10c1738', \
'1.0.4' : '23867540ec32e63ce71d6ecc105278d2', \
}
# extra help message

2
src/.gitignore vendored
View File

@ -1801,4 +1801,6 @@
/pair_smtbq.h
/pair_vashishta*.cpp
/pair_vashishta*.h
/pair_lj_pirani.cpp
/pair_lj_pirani.h

7
src/BPM/bond_bpm.cpp
View File

@ -397,6 +397,7 @@ double BondBPM::equilibrium_distance(int /*i*/)
void BondBPM::write_restart(FILE *fp)
{
fwrite(&overlay_flag, sizeof(int), 1, fp);
fwrite(&break_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
@ -405,8 +406,12 @@ void BondBPM::write_restart(FILE *fp)
void BondBPM::read_restart(FILE *fp)
{
if (comm->me == 0) utils::sfread(FLERR, &overlay_flag, sizeof(int), 1, fp, nullptr, error);
if (comm->me == 0) {
utils::sfread(FLERR, &overlay_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &break_flag, sizeof(int), 1, fp, nullptr, error);
}
MPI_Bcast(&overlay_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&break_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */

11
src/BPM/bond_bpm_rotational.cpp
View File

@ -28,6 +28,7 @@
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
@ -483,6 +484,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
int newton_bond = force->newton_bond;
double **bondstore = fix_bond_history->bondstore;
const bool allow_breaks = (update->setupflag == 0) && break_flag;
for (n = 0; n < nbondlist; n++) {
@ -527,7 +529,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
breaking = elastic_forces(i1, i2, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
torque1on2, torque2on1);
if ((breaking >= 1.0) && break_flag) {
if ((breaking >= 1.0) && allow_breaks) {
bondlist[n][2] = 0;
process_broken(i1, i2);
continue;
@ -763,6 +765,7 @@ void BondBPMRotational::read_restart(FILE *fp)
void BondBPMRotational::write_restart_settings(FILE *fp)
{
fwrite(&smooth_flag, sizeof(int), 1, fp);
fwrite(&normalize_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
@ -771,8 +774,12 @@ void BondBPMRotational::write_restart_settings(FILE *fp)
void BondBPMRotational::read_restart_settings(FILE *fp)
{
if (comm->me == 0) utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error);
if (comm->me == 0) {
utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &normalize_flag, sizeof(int), 1, fp, nullptr, error);
}
MPI_Bcast(&smooth_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&normalize_flag, 1, MPI_INT, 0, world);
}
/* ---------------------------------------------------------------------- */

4
src/BPM/bond_bpm_spring.cpp
View File

@ -26,6 +26,7 @@
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
@ -218,6 +219,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
double invdim = 1.0 / dim;
double **bondstore = fix_bond_history->bondstore;
const bool allow_breaks = (update->setupflag == 0) && break_flag;
for (n = 0; n < nbondlist; n++) {
@ -249,7 +251,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
r = sqrt(rsq);
e = (r - r0) / r0;
if ((fabs(e) > ecrit[type]) && break_flag) {
if ((fabs(e) > ecrit[type]) && allow_breaks) {
bondlist[n][2] = 0;
process_broken(i1, i2);

4
src/BPM/bond_bpm_spring_plastic.cpp
View File

@ -26,6 +26,7 @@
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
@ -185,6 +186,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag)
int newton_bond = force->newton_bond;
double **bondstore = fix_bond_history->bondstore;
const bool allow_breaks = (update->setupflag == 0) && break_flag;
for (n = 0; n < nbondlist; n++) {
@ -217,7 +219,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag)
r = sqrt(rsq);
e = (r - r0) / r0;
if ((fabs(e) > ecrit[type]) && break_flag) {
if ((fabs(e) > ecrit[type]) && allow_breaks) {
bondlist[n][2] = 0;
process_broken(i1, i2);
continue;

31
src/BPM/fix_update_special_bonds.cpp
View File

@ -26,6 +26,7 @@
#include "neigh_list.h"
#include "neighbor.h"
#include <set>
#include <utility>
using namespace LAMMPS_NS;
@ -267,18 +268,40 @@ void FixUpdateSpecialBonds::post_run()
void FixUpdateSpecialBonds::add_broken_bond(int i, int j)
{
auto tag_pair = std::make_pair(atom->tag[i], atom->tag[j]);
tagint *tag = atom->tag;
int mintag = MIN(tag[i], tag[j]);
int maxtag = MAX(tag[i], tag[j]);
auto tag_pair = std::make_pair(mintag, maxtag);
new_broken_pairs.push_back(tag_pair);
broken_pairs.push_back(tag_pair);
// cancel out if created->destroyed before nlist rebuild
// however, still cannot break + create in the same timestep
auto const &it = created_pairs.find(tag_pair);
if (it != created_pairs.end()) {
created_pairs.erase(it);
} else {
broken_pairs.insert(tag_pair);
}
}
/* ---------------------------------------------------------------------- */
void FixUpdateSpecialBonds::add_created_bond(int i, int j)
{
auto tag_pair = std::make_pair(atom->tag[i], atom->tag[j]);
tagint *tag = atom->tag;
int mintag = MIN(tag[i], tag[j]);
int maxtag = MAX(tag[i], tag[j]);
auto tag_pair = std::make_pair(mintag, maxtag);
new_created_pairs.push_back(tag_pair);
created_pairs.push_back(tag_pair);
// cancel out if destroyed->created before nlist rebuild
// however, still cannot break + create in the same timestep
auto const &it = broken_pairs.find(tag_pair);
if (it != broken_pairs.end()) {
broken_pairs.erase(it);
} else {
created_pairs.insert(tag_pair);
}
}
/* ----------------------------------------------------------------------

13
src/BPM/fix_update_special_bonds.h
View File

@ -22,6 +22,7 @@ FixStyle(UPDATE_SPECIAL_BONDS,FixUpdateSpecialBonds);
#include "fix.h"
#include <set>
#include <utility>
namespace LAMMPS_NS {
@ -39,15 +40,15 @@ class FixUpdateSpecialBonds : public Fix {
void write_restart(FILE *) override;
protected:
// Create two arrays to store bonds broken this timestep (new)
// and since the last neighbor list build
// Create array to store bonds broken this timestep (new)
// and a set for those broken since the last neighbor list build
std::vector<std::pair<tagint, tagint>> new_broken_pairs;
std::vector<std::pair<tagint, tagint>> broken_pairs;
std::set<std::pair<tagint, tagint>> broken_pairs;
// Create two arrays to store newly created this timestep (new)
// and since the last neighbor list build
// Create arrays to store newly created this timestep (new)
// and a set for those created since the last neighbor list build
std::vector<std::pair<tagint, tagint>> new_created_pairs;
std::vector<std::pair<tagint, tagint>> created_pairs;
std::set<std::pair<tagint, tagint>> created_pairs;
};
} // namespace LAMMPS_NS

886
src/EXTRA-PAIR/pair_lj_pirani.cpp Normal file
View File

@ -0,0 +1,886 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mateo Rodríguez (mateorsuarez@gmail.com) (IFF-CSIC)
Work done at the Molecular Interactions Group (INTERMOL) of the
Fundamental Physics Institute (http://intermol.iff.csic.es/).
Optimization of the code: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "pair_lj_pirani.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_special.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using MathSpecial::square;
/* ---------------------------------------------------------------------- */
PairLJPirani::PairLJPirani(LAMMPS *lmp) : Pair(lmp), cut_respa(nullptr)
{
respa_enable = 1;
born_matrix_enable = 0;
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairLJPirani::~PairLJPirani()
{
if (copymode) return;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(alpha);
memory->destroy(beta);
memory->destroy(gamma);
memory->destroy(rm);
memory->destroy(epsilon);
memory->destroy(offset);
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJPirani::allocate()
{
allocated = 1;
int n = atom->ntypes + 1;
memory->create(setflag, n, n, "pair:setflag");
for (int i = 1; i < n; i++)
for (int j = i; j < n; j++) setflag[i][j] = 0;
memory->create(cutsq, n, n, "pair:cutsq");
memory->create(cut, n, n, "pair:cut");
memory->create(alpha, n, n, "pair:alpha");
memory->create(beta, n, n, "pair:beta");
memory->create(gamma, n, n, "pair:gamma");
memory->create(rm, n, n, "pair:rm");
memory->create(epsilon, n, n, "pair:epsilon");
memory->create(offset, n, n, "pair:offset");
}
/* ---------------------------------------------------------------------- */
void PairLJPirani::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, rx, n_x;
double pow_rx_n_x, pow_rx_gamma;
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
double ilj1, ilj2;
double fxtmp, fytmp, fztmp;
evdwl = 0.0;
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
fxtmp = fytmp = fztmp = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
rx = r / rm[itype][jtype];
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
pow_rx_n_x = pow(1.0 / rx, n_x);
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
// F = -dV/dr
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
fpair = factor_lj * forceilj / r; // F_x = -x/r * dV/dr (chain rule)
fxtmp += delx * fpair;
fytmp += dely * fpair;
fztmp += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (eflag) {
ilj1 = epsilon[itype][jtype] * gamma[itype][jtype] * pow(1 / rx, n_x) /
(n_x - gamma[itype][jtype]);
ilj2 = -epsilon[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) /
(n_x - gamma[itype][jtype]);
evdwl = ilj1 + ilj2 - offset[itype][jtype];
evdwl *= factor_lj;
}
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void PairLJPirani::compute_inner()
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
double rsq, factor_lj, rsw;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, rx, n_x;
double pow_rx_n_x, pow_rx_gamma;
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
double fxtmp, fytmp, fztmp;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
inum = list->inum_inner;
ilist = list->ilist_inner;
numneigh = list->numneigh_inner;
firstneigh = list->firstneigh_inner;
double cut_out_on = cut_respa[0];
double cut_out_off = cut_respa[1];
double cut_out_diff = cut_out_off - cut_out_on;
double cut_out_on_sq = cut_out_on * cut_out_on;
double cut_out_off_sq = cut_out_off * cut_out_off;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
fxtmp = fytmp = fztmp = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cut_out_off_sq) {
jtype = type[j];
r = sqrt(rsq);
rx = r / rm[itype][jtype];
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
pow_rx_n_x = pow(1.0 / rx, n_x);
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
// F = -dV/dr
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
fpair = factor_lj * forceilj / r; // F_x = -x/r * dV/dr (chain rule)
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
fpair *= 1.0 - rsw * rsw * (3.0 - 2.0 * rsw);
}
fxtmp += delx * fpair;
fytmp += dely * fpair;
fztmp += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
}
/* ---------------------------------------------------------------------- */
void PairLJPirani::compute_middle()
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
double rsq, factor_lj, rsw;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, rx, n_x;
double pow_rx_n_x, pow_rx_gamma;
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
double fxtmp, fytmp, fztmp;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
inum = list->inum_middle;
ilist = list->ilist_middle;
numneigh = list->numneigh_middle;
firstneigh = list->firstneigh_middle;
double cut_in_off = cut_respa[0];
double cut_in_on = cut_respa[1];
double cut_out_on = cut_respa[2];
double cut_out_off = cut_respa[3];
double cut_in_diff = cut_in_on - cut_in_off;
double cut_out_diff = cut_out_off - cut_out_on;
double cut_in_off_sq = cut_in_off * cut_in_off;
double cut_in_on_sq = cut_in_on * cut_in_on;
double cut_out_on_sq = cut_out_on * cut_out_on;
double cut_out_off_sq = cut_out_off * cut_out_off;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
fxtmp = fytmp = fztmp = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
jtype = type[j];
r = sqrt(rsq);
rx = r / rm[itype][jtype];
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
pow_rx_n_x = pow(1.0 / rx, n_x);
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
// F = -dV/dr
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
fpair = factor_lj * forceilj / r;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
}
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0);
}
fxtmp += delx * fpair;
fytmp += dely * fpair;
fztmp += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
}
/* ---------------------------------------------------------------------- */
void PairLJPirani::compute_outer(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq, factor_lj, rsw;
int *ilist, *jlist, *numneigh, **firstneigh;
double r, rx, n_x;
double pow_rx_n_x, pow_rx_gamma;
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
double ilj1, ilj2;
double fxtmp, fytmp, fztmp;
evdwl = 0.0;
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double cut_in_off = cut_respa[2];
double cut_in_on = cut_respa[3];
double cut_in_diff = cut_in_on - cut_in_off;
double cut_in_off_sq = cut_in_off * cut_in_off;
double cut_in_on_sq = cut_in_on * cut_in_on;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
fxtmp = fytmp = fztmp = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
if (rsq > cut_in_off_sq) {
r = sqrt(rsq);
rx = r / rm[itype][jtype];
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
pow_rx_n_x = pow(1.0 / rx, n_x);
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
filj6 =
+1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
// F = -dV/dr
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
fpair = factor_lj * forceilj / r;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
}
fxtmp += delx * fpair;
fytmp += dely * fpair;
fztmp += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
}
if (eflag) {
r = sqrt(rsq);
rx = r / rm[itype][jtype];
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
ilj1 = epsilon[itype][jtype] * gamma[itype][jtype] * pow(1 / rx, n_x) /
(n_x - gamma[itype][jtype]);
ilj2 = -epsilon[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) /
(n_x - gamma[itype][jtype]);
evdwl = ilj1 + ilj2 - offset[itype][jtype];
evdwl *= factor_lj;
}
if (vflag) {
if (rsq <= cut_in_off_sq) {
r = sqrt(rsq);
rx = r / rm[itype][jtype];
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
pow_rx_n_x = pow(1.0 / rx, n_x);
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
filj5 =
-1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
filj6 =
+1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
// F = -dV/dr
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
fpair = factor_lj * forceilj / r;
} else if (rsq < cut_in_on_sq)
fpair = factor_lj * forceilj / r;
}
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
}
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJPirani::settings(int narg, char **arg)
{
if (narg != 1)
error->all(FLERR, "Pair style ilj/cut must have exactly one argument: cutoff distance");
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
// reset cutoffs that have been explicitly set
if (allocated) {
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global;
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
/*
7 or 8 coefficients: 5 for the ILJ, 2 for the pair, 1 for the cutoff (optional)
*/
void PairLJPirani::coeff(int narg, char **arg)
{
if (narg < 7 || narg > 8) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double alpha_one = utils::numeric(FLERR, arg[2], false, lmp);
double beta_one = utils::numeric(FLERR, arg[3], false, lmp);
double gamma_one = utils::numeric(FLERR, arg[4], false, lmp);
double rm_one = utils::numeric(FLERR, arg[5], false, lmp);
double epsilon_one = utils::numeric(FLERR, arg[6], false, lmp);
double cut_one = cut_global;
if (narg == 8) cut_one = utils::numeric(FLERR, arg[7], false, lmp);
if (rm_one <= 0.0 || epsilon_one < 0.0 || gamma_one <= 0.0)
error->all(FLERR, "Illegal ILJ coefficients");
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
alpha[i][j] = alpha_one;
beta[i][j] = beta_one;
gamma[i][j] = gamma_one;
rm[i][j] = rm_one;
epsilon[i][j] = epsilon_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
// Initialize symmetric entries
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
alpha[j][i] = alpha[i][j];
beta[j][i] = beta[i][j];
gamma[j][i] = gamma[i][j];
rm[j][i] = rm[i][j];
epsilon[j][i] = epsilon[i][j];
cut[j][i] = cut[i][j];
setflag[j][i] = setflag[i][j];
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJPirani::init_style()
{
// request regular or rRESPA neighbor list
int list_style = NeighConst::REQ_DEFAULT;
if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) {
auto respa = dynamic_cast<Respa *>(update->integrate);
if (respa->level_inner >= 0) list_style = NeighConst::REQ_RESPA_INOUT;
if (respa->level_middle >= 0) list_style = NeighConst::REQ_RESPA_ALL;
}
neighbor->add_request(this, list_style);
// set rRESPA cutoffs
if (utils::strmatch(update->integrate_style, "^respa") &&
(dynamic_cast<Respa *>(update->integrate))->level_inner >= 0)
cut_respa = (dynamic_cast<Respa *>(update->integrate))->cutoff;
else
cut_respa = nullptr;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJPirani::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
if (offset_flag && (cut[i][j] > 0.0)) {
double r = cut[i][j] / rm[i][j];
double nx = alpha[i][j] * r * r + beta[i][j];
offset[i][j] = epsilon[i][j] *
((gamma[i][j] / (nx - gamma[i][j])) * pow(1 / r, nx) -
(nx / (nx - gamma[i][j])) * pow(1 / r, gamma[i][j]));
} else
offset[i][j] = 0.0;
alpha[j][i] = alpha[i][j];
beta[j][i] = beta[i][j];
gamma[j][i] = gamma[i][j];
rm[j][i] = rm[i][j];
epsilon[j][i] = epsilon[i][j];
offset[j][i] = offset[i][j];
// check interior rRESPA cutoff
if (cut_respa && cut[i][j] < cut_respa[3])
error->all(FLERR, "Pair cutoff < Respa interior cutoff");
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJPirani::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&alpha[i][j], sizeof(double), 1, fp);
fwrite(&beta[i][j], sizeof(double), 1, fp);
fwrite(&gamma[i][j], sizeof(double), 1, fp);
fwrite(&rm[i][j], sizeof(double), 1, fp);
fwrite(&epsilon[i][j], sizeof(double), 1, fp);
fwrite(&cut[i][j], sizeof(double), 1, fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJPirani::write_restart_settings(FILE *fp)
{
fwrite(&cut_global, sizeof(double), 1, fp);
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJPirani::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR, &alpha[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &beta[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &gamma[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &rm[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&alpha[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&beta[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&gamma[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&rm[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJPirani::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
}
MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairLJPirani::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp, "%d %g %g %g %g %g\n", i, alpha[i][i], beta[i][i], gamma[i][i], rm[i][i],
epsilon[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairLJPirani::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp, "%d %d %g %g %g %g %g %g\n", i, j, alpha[i][j], beta[i][j], gamma[i][j], rm[i][j],
epsilon[i][j], cut[i][j]);
}
/* ---------------------------------------------------------------------- */
double PairLJPirani::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double /*factor_coul*/, double factor_lj, double &fforce)
{
double r, rx, n_x, filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
double ilj1, ilj2;
r = sqrt(rsq);
rx = r / rm[itype][jtype];
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow(1 / rx, n_x) /
(pow(n_x - gamma[itype][jtype], 2.0) * rm[itype][jtype]);
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow(1 / rx, gamma[itype][jtype]) /
(pow(n_x - gamma[itype][jtype], 2.0) * rm[itype][jtype]);
filj3 = -2.0 * alpha[itype][jtype] * rx * pow(1 / rx, gamma[itype][jtype]) /
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) * log(1 / rx) *
pow(1 / rx, n_x) / (n_x - gamma[itype][jtype]);
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow(1 / rx, n_x) / (r * (n_x - gamma[itype][jtype]));
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) /
(r * (n_x - gamma[itype][jtype]));
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
fforce = factor_lj * forceilj / r;
ilj1 =
epsilon[itype][jtype] * gamma[itype][jtype] * pow(1 / rx, n_x) / (n_x - gamma[itype][jtype]);
ilj2 =
-epsilon[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) / (n_x - gamma[itype][jtype]);
return factor_lj * (ilj1 + ilj2 - offset[itype][jtype]);
}
/* ---------------------------------------------------------------------- */
void *PairLJPirani::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str, "alpha") == 0) return (void *) alpha;
if (strcmp(str, "beta") == 0) return (void *) beta;
if (strcmp(str, "gamma") == 0) return (void *) gamma;
if (strcmp(str, "rm") == 0) return (void *) rm;
if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
return nullptr;
}

60
src/EXTRA-PAIR/pair_lj_pirani.h Normal file
View File

@ -0,0 +1,60 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(lj/pirani,PairLJPirani);
// clang-format on
#else
#ifndef LMP_PAIR_LJ_PIRANI_H
#define LMP_PAIR_LJ_PIRANI_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJPirani : public Pair {
public:
PairLJPirani(class LAMMPS *);
virtual ~PairLJPirani() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
void write_data(FILE *) override;
void write_data_all(FILE *) override;
double single(int, int, int, int, double, double, double, double &) override;
void *extract(const char *, int &) override;
void compute_inner() override;
void compute_middle() override;
void compute_outer(int, int) override;
protected:
double cut_global;
double **cut;
double **alpha, **beta, **gamma, **rm, **epsilon;
double **offset;
double *cut_respa;
virtual void allocate();
};
} // namespace LAMMPS_NS
#endif
#endif

65
src/MANYBODY/pair_eam.cpp
View File

@ -39,6 +39,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
{
restartinfo = 0;
manybody_flag = 1;
atomic_energy_enable = 1;
embedstep = -1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
@ -340,6 +341,70 @@ void PairEAM::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
}
/*********************************************************************
* Calculates the atomic energy of atom i
*********************************************************************/
double PairEAM::compute_atomic_energy(int i, NeighList *neighborList)
{
double p;
int m;
double* coeff;
double Ei = 0.0;
double rhoi = 0.0;
double xi = atom->x[i][0];
double yi = atom->x[i][1];
double zi = atom->x[i][2];
int itype = atom->type[i];
// loop over all neighbors of the selected atom.
int* jlist = neighborList->firstneigh[i];
int jnum = neighborList->numneigh[i];
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
double delx = xi - atom->x[j][0];
double dely = yi - atom->x[j][1];
double delz = zi - atom->x[j][2];
double rsq = delx*delx + dely*dely + delz*delz;
if(rsq >= cutforcesq) continue;
int jtype = atom->type[j];
double r = sqrt(rsq);
p = r * rdr + 1.0;
m = static_cast<int>(p);
m = MIN(m, nr - 1);
p -= m;
p = MIN(p, 1.0);
// sum pair energy ij
// divide by 2 to avoid double counting energy
coeff = z2r_spline[type2z2r[jtype][itype]][m];
double z2 = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
Ei += 0.5*z2 / r;
// sum rho_ij to rho_i
coeff = rhor_spline[type2rhor[jtype][itype]][m];
rhoi += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
}
// compute the change in embedding energy of atom i.
p = rhoi * rdrho + 1.0;
m = static_cast<int>(p);
m = MAX(1, MIN(m, nrho - 1));
p -= m;
p = MIN(p, 1.0);
coeff = frho_spline[type2frho[itype]][m];
Ei += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
return Ei;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */

1
src/MANYBODY/pair_eam.h
View File

@ -47,6 +47,7 @@ class PairEAM : public Pair {
PairEAM(class LAMMPS *);
~PairEAM() override;
void compute(int, int) override;
double compute_atomic_energy(int, NeighList *) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;

1
src/MANYBODY/pair_extep.cpp
View File

@ -161,7 +161,6 @@ void PairExTeP::SR_neigh()
}
}
}
//printf("SR_neigh : N[%d] = %f\n",i,N[i]);
ipage->vgot(n);
if (ipage->status())

2
src/MANYBODY/pair_gw.h
View File

@ -94,8 +94,6 @@ class PairGW : public Pair {
const double gw_d = param->d * param->d;
const double hcth = param->h - costheta;
//printf("gw_gijk: gw_c=%f gw_d=%f hcth=%f=%f-%f\n", gw_c, gw_d, hcth, param->h, costheta);
return param->gamma * (1.0 + gw_c / gw_d - gw_c / (gw_d + hcth * hcth));
}

5
src/MC/fix_bond_break.cpp
View File

@ -439,11 +439,6 @@ void FixBondBreak::update_topology()
ndihedrals = 0;
nimpropers = 0;
//printf("NBREAK %d: ",nbreak);
//for (i = 0; i < nbreak; i++)
// printf(" %d %d,",broken[i][0],broken[i][1]);
//printf("\n");
for (i = 0; i < nlocal; i++) {
influenced = 0;
slist = special[i];

137
src/MC/fix_sgcmc.cpp
View File

@ -124,6 +124,7 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char **
// Default values for optional parameters (where applicable).
numSamplingWindowMoves = 8;
seed = 324234;
atomicenergyflag = 0;
// Parse extra/optional parameters
while (iarg < narg) {
@ -175,7 +176,7 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char **
targetConcentration[i] = utils::numeric(FLERR, arg[iarg], false, lmp);
targetConcentration[1] -= targetConcentration[i];
}
for (int i = 1; i <= atom->ntypes; i++, iarg++) {
for (int i = 1; i <= atom->ntypes; i++) {
if ((targetConcentration[i] < 0.0) || (targetConcentration[i] > 1.0))
error->all(FLERR, "Target concentration {} for species {} is out of range",
targetConcentration[i], i);
@ -183,7 +184,6 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char **
utils::logmesg(lmp, " SGC - Target concentration of species {}: {}\n", i,
targetConcentration[i]);
}
} else if (strcmp(arg[iarg], "serial") == 0) {
// Switch off second rejection.
serialMode = true;
@ -193,6 +193,12 @@ FixSemiGrandCanonicalMC::FixSemiGrandCanonicalMC(LAMMPS *_lmp, int narg, char **
if (comm->nprocs != 1)
error->all(FLERR, "Cannot use serial mode Monte Carlo in a parallel simulation.");
} else if (strcmp(arg[iarg],"atomic/energy") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, fmt::format("fix {} atomic/energy", style), error);
atomicenergyflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else {
error->all(FLERR, "Unknown fix sgcmc keyword: {}", arg[iarg]);
}
@ -243,25 +249,40 @@ void FixSemiGrandCanonicalMC::init()
if (modify->get_fix_by_style("sgcmc").size() > 1)
error->all(FLERR, "More than one fix sgcmc defined.");
// Save a pointer to the EAM potential.
pairEAM = dynamic_cast<PairEAM*>(force->pair);
if (!pairEAM) {
if (comm->me == 0)
utils::logmesg(lmp, " SGC - Using naive total energy calculation for MC -> SLOW!\n");
if (atomicenergyflag) {
// Save a pointer to the EAM potential.
if (comm->me == 0) utils::logmesg(lmp, " SGC - Using atomic energy method for SGCMC\n");
if (!force->pair->atomic_energy_enable) {
error->all(FLERR, "SGC - Pair style does not support atomic energy method");
}
} else {
// Save a pointer to the EAM potential.
pairEAM = dynamic_cast<PairEAM*>(force->pair);
if (!pairEAM) {
if (comm->nprocs > 1)
error->all(FLERR, "Can not run fix sgcmc with naive total energy calculation "
"and more than one MPI process.");
if (comm->me == 0)
utils::logmesg(lmp, " SGC - Using naive total energy calculation for MC -> SLOW!\n");
// Get reference to a compute that will provide the total energy of the system.
// This is needed by computeTotalEnergy().
compute_pe = modify->get_compute_by_id("thermo_pe");
if (comm->nprocs > 1)
error->all(FLERR, "Can not run fix sgcmc with naive total energy calculation "
"and more than one MPI process.");
// Get reference to a compute that will provide the total energy of the system.
// This is needed by computeTotalEnergy().
compute_pe = modify->get_compute_by_id("thermo_pe");
}
}
interactionRadius = force->pair->cutforce;
if (comm->me == 0) utils::logmesg(lmp, " SGC - Interaction radius: {}\n", interactionRadius);
// This fix needs a full neighbor list.
neighbor->add_request(this, NeighConst::REQ_FULL);
if (atomicenergyflag)
// for atomic energy method, need ghost neighbors
neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_GHOST);
else
neighbor->add_request(this, NeighConst::REQ_FULL);
// Count local number of atoms from each species.
const int *type = atom->type;
@ -379,11 +400,18 @@ void FixSemiGrandCanonicalMC::doMC()
deltaN[newSpecies] = +1;
// Compute the energy difference that swapping this atom would cost or gain.
if (pairEAM) {
deltaE = computeEnergyChangeEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies);
// Atomic energy method:
if(atomicenergyflag) {
deltaE = computeEnergyChangeEatom(selectedAtom, oldSpecies, newSpecies);
} else {
// Generic case:
deltaE = computeEnergyChangeGeneric(selectedAtom, oldSpecies, newSpecies);
// EAM method:
if (pairEAM) {
deltaE = computeEnergyChangeEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies);
} else {
// Generic method
deltaE = computeEnergyChangeGeneric(selectedAtom, oldSpecies, newSpecies);
}
}
// Perform inner MC acceptance test.
@ -440,10 +468,14 @@ void FixSemiGrandCanonicalMC::doMC()
// Make accepted atom swap permanent.
if (selectedAtom >= 0) {
if (pairEAM)
flipAtomEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies);
else
flipAtomGeneric(selectedAtom, oldSpecies, newSpecies);
if(atomicenergyflag) {
flipAtomEatom(selectedAtom, oldSpecies, newSpecies);
} else {
if (pairEAM)
flipAtomEAM(selectedAtom, selectedAtomNL, oldSpecies, newSpecies);
else
flipAtomGeneric(selectedAtom, oldSpecies, newSpecies);
}
nAcceptedSwapsLocal++;
} else {
nRejectedSwapsLocal++;
@ -989,6 +1021,67 @@ void FixSemiGrandCanonicalMC::flipAtomGeneric(int flipAtom, int oldSpecies, int
changedAtoms[flipAtom] = true;
}
/*********************************************************************
* Calculates the change in energy that swapping the given
* atom would produce. This routine uses a per-atom energy calculation
*********************************************************************/
double FixSemiGrandCanonicalMC::computeEnergyChangeEatom(int flipAtom, int oldSpecies, int newSpecies)
{
double Eold, Enew, deltaE;
// Calculate old atomic energy of selected atom
Eold = force->pair->compute_atomic_energy(flipAtom, neighborList);
// calculate the old per-atom energy of neighbors
int* jlist = neighborList->firstneigh[flipAtom];
int jnum = neighborList->numneigh[flipAtom];
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
Eold += force->pair->compute_atomic_energy(j, neighborList);
}
// Calculate new per-atom energy of selected atom
atom->type[flipAtom] = newSpecies;
Enew = force->pair->compute_atomic_energy(flipAtom, neighborList);
// calculate the new per-atom energy of neighbors
for(int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
Enew += force->pair->compute_atomic_energy(j, neighborList);
}
atom->type[flipAtom] = oldSpecies;
deltaE = Enew - Eold;
return deltaE;
}
/*********************************************************************
* Flips the type of one atom.
* This routine is for the per-atom energy case.
*********************************************************************/
void FixSemiGrandCanonicalMC::flipAtomEatom(int flipAtom, int oldSpecies, int newSpecies)
{
atom->type[flipAtom] = newSpecies;
// Rescale particle velocity vector to conserve kinetic energy.
double vScaleFactor = sqrt(atom->mass[oldSpecies] / atom->mass[newSpecies]);
atom->v[flipAtom][0] *= vScaleFactor;
atom->v[flipAtom][1] *= vScaleFactor;
atom->v[flipAtom][2] *= vScaleFactor;
changedAtoms[flipAtom] = true;
}
/*********************************************************************
* Lets the fix report one of its internal state variables to LAMMPS.
*********************************************************************/

11
src/MC/fix_sgcmc.h
View File

@ -75,6 +75,10 @@ class FixSemiGrandCanonicalMC : public Fix {
// This routine should only be used for debugging purposes.
double computeEnergyChangeGeneric(int flipAtom, int oldSpecies, int newSpecies);
// Calculates the change in energy that swapping the given atom would produce.
// This uses the atomic energy method
double computeEnergyChangeEatom(int flipAtom, int oldSpecies, int newSpecies);
// Lets LAMMPS calculate the total potential energy of the system.
double computeTotalEnergy();
@ -86,6 +90,10 @@ class FixSemiGrandCanonicalMC : public Fix {
// This routine is for the generic case.
void flipAtomGeneric(int flipAtom, int oldSpecies, int newSpecies);
// Flips the type of one atom.
// This routine is for the atomic energy method
void flipAtomEatom(int flipAtom, int oldSpecies, int newSpecies);
// Transfers the locally changed electron densities and atom types to the neighbors.
void communicateRhoAndTypes();
@ -181,6 +189,9 @@ class FixSemiGrandCanonicalMC : public Fix {
// A compute used to compute the total potential energy of the system.
class Compute *compute_pe;
// Indicate whether or not atomic energy method is used
int atomicenergyflag;
};
} // namespace LAMMPS_NS

194
src/OPENMP/pair_lj_pirani_omp.cpp Normal file
View File

@ -0,0 +1,194 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "pair_lj_pirani_omp.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "math_special.h"
#include "neigh_list.h"
#include "suffix.h"
#include "omp_compat.h"
using namespace LAMMPS_NS;
using MathSpecial::square;
/* ---------------------------------------------------------------------- */
PairLJPiraniOMP::PairLJPiraniOMP(LAMMPS *lmp) : PairLJPirani(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
respa_enable = 0;
cut_respa = nullptr;
}
/* ---------------------------------------------------------------------- */
void PairLJPiraniOMP::compute(int eflag, int vflag)
{
ev_init(eflag, vflag);
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
const int inum = list->inum;
#if defined(_OPENMP)
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag)
#endif
{
int ifrom, ito, tid;
loop_setup_thr(ifrom, ito, tid, inum, nthreads);
ThrData *thr = fix->get_thr(tid);
thr->timer(Timer::START);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
if (evflag) {
if (eflag) {
if (force->newton_pair)
eval<1, 1, 1>(ifrom, ito, thr);
else
eval<1, 1, 0>(ifrom, ito, thr);
} else {
if (force->newton_pair)
eval<1, 0, 1>(ifrom, ito, thr);
else
eval<1, 0, 0>(ifrom, ito, thr);
}
} else {
if (force->newton_pair)
eval<0, 0, 1>(ifrom, ito, thr);
else
eval<0, 0, 0>(ifrom, ito, thr);
}
thr->timer(Timer::PAIR);
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
}
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void PairLJPiraniOMP::eval(int iifrom, int iito, ThrData *const thr)
{
const auto *_noalias const x = (dbl3_t *) atom->x[0];
auto *_noalias const f = (dbl3_t *) thr->get_f()[0];
const int *_noalias const type = atom->type;
const double *_noalias const special_lj = force->special_lj;
const int *_noalias const ilist = list->ilist;
const int *_noalias const numneigh = list->numneigh;
const int *const *const firstneigh = list->firstneigh;
double xtmp, ytmp, ztmp, delx, dely, delz, fxtmp, fytmp, fztmp;
double rsq, forceilj, factor_lj, evdwl, fpair;
const int nlocal = atom->nlocal;
int j, jj, jnum, jtype;
evdwl = 0.0;
// loop over neighbors of my atoms
for (int ii = iifrom; ii < iito; ++ii) {
const int i = ilist[ii];
const int itype = type[i];
const int *_noalias const jlist = firstneigh[i];
const double *_noalias const cutsqi = cutsq[itype];
const double *_noalias const offseti = offset[itype];
const double *_noalias const alphai = alpha[itype];
const double *_noalias const betai = beta[itype];
const double *_noalias const gammai = gamma[itype];
const double *_noalias const rmi = rm[itype];
const double *_noalias const epsiloni = epsilon[itype];
xtmp = x[i].x;
ytmp = x[i].y;
ztmp = x[i].z;
jnum = numneigh[i];
fxtmp = fytmp = fztmp = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j].x;
dely = ytmp - x[j].y;
delz = ztmp - x[j].z;
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsqi[jtype]) {
const double r = sqrt(rsq);
const double rx = r / rmi[jtype];
const double n_x = alphai[jtype] * rx * rx + betai[jtype];
const double pow_rx_n_x = pow(1.0 / rx, n_x);
const double pow_rx_gamma = pow(1.0 / rx, gammai[jtype]);
double filj1 = -2.0 * alphai[jtype] * gammai[jtype] * rx * pow_rx_n_x /
(square(n_x - gammai[jtype]) * rmi[jtype]);
double filj2 = 2.0 * alphai[jtype] * rx * n_x * pow_rx_gamma /
(square(n_x - gammai[jtype]) * rmi[jtype]);
double filj3 =
-2.0 * alphai[jtype] * rx * pow_rx_gamma / (rmi[jtype] * (n_x - gammai[jtype]));
double filj4 = 2.0 * alphai[jtype] * gammai[jtype] * (rx / rmi[jtype]) * log(1 / rx) *
pow_rx_n_x / (n_x - gammai[jtype]);
double filj5 = -1.0 * gammai[jtype] * n_x * pow_rx_n_x / (r * (n_x - gammai[jtype]));
double filj6 = 1.0 * gammai[jtype] * n_x * pow_rx_gamma / (r * (n_x - gammai[jtype]));
// F = -dV/dr
forceilj = -epsiloni[jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
fpair = factor_lj * forceilj / r; // F_x = -x/r * dV/dr (chain rule)
fxtmp += delx * fpair;
fytmp += dely * fpair;
fztmp += delz * fpair;
if (NEWTON_PAIR || j < nlocal) {
f[j].x -= delx * fpair;
f[j].y -= dely * fpair;
f[j].z -= delz * fpair;
}
if (EFLAG) {
double ilj1 = epsiloni[jtype] * gammai[jtype] * pow(1 / rx, n_x) / (n_x - gammai[jtype]);
double ilj2 = -epsiloni[jtype] * n_x * pow(1 / rx, gammai[jtype]) / (n_x - gammai[jtype]);
evdwl = ilj1 + ilj2 - offseti[jtype];
evdwl *= factor_lj;
}
if (EVFLAG)
ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz, thr);
}
}
f[i].x += fxtmp;
f[i].y += fytmp;
f[i].z += fztmp;
}
}
/* ---------------------------------------------------------------------- */
double PairLJPiraniOMP::memory_usage()
{
double bytes = memory_usage_thr();
bytes += PairLJPirani::memory_usage();
return bytes;
}

48
src/OPENMP/pair_lj_pirani_omp.h Normal file
View File

@ -0,0 +1,48 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(lj/pirani/omp,PairLJPiraniOMP);
// clang-format on
#else
#ifndef LMP_PAIR_LJ_PIRANI_OMP_H
#define LMP_PAIR_LJ_PIRANI_OMP_H
#include "pair_lj_pirani.h"
#include "thr_omp.h"
namespace LAMMPS_NS {
class PairLJPiraniOMP : public PairLJPirani, public ThrOMP {
public:
PairLJPiraniOMP(class LAMMPS *);
void compute(int, int) override;
double memory_usage() override;
private:
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
void eval(int ifrom, int ito, ThrData *const thr);
};
} // namespace LAMMPS_NS
#endif
#endif

2
src/REAXFF/reaxff_tool_box.cpp
View File

@ -101,6 +101,6 @@ void sfree(void *ptr)
if (ptr == nullptr) return;
free(ptr);
ptr = nullptr;
ptr = nullptr; // NOTE: this has no effect
}
} // namespace ReaxFF

4
src/RHEO/bond_rheo_shell.cpp
View File

@ -44,7 +44,7 @@ using namespace RHEO_NS;
/* ---------------------------------------------------------------------- */
BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr),
BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr),
id_fix(nullptr), compute_surface(nullptr)
{
partial_flag = 1;
@ -71,7 +71,6 @@ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
modify->add_fix(fmt::format("{} all property/atom i_shell_nbond", id_fix));
index_nb = atom->find_custom("shell_nbond", tmp1, tmp2);
}
nbond = atom->ivector[index_nb];
//Store non-persistent per atom quantities, intermediate
@ -181,6 +180,7 @@ void BondRHEOShell::compute(int eflag, int vflag)
double **v = atom->v;
double **f = atom->f;
tagint *tag = atom->tag;
int *nbond = atom->ivector[index_nb];
int *status = atom->rheo_status;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;

2
src/RHEO/bond_rheo_shell.h
View File

@ -45,7 +45,7 @@ class BondRHEOShell : public BondBPM {
double *k, *ecrit, *gamma;
double tform, rmax, rsurf;
int *dbond, *nbond;
int *dbond;
int index_nb, nmax_store;
char *id_fix;

24
src/RHEO/compute_rheo_interface.cpp
View File

@ -39,12 +39,10 @@ using namespace LAMMPS_NS;
using namespace RHEO_NS;
using namespace MathExtra;
static constexpr double EPSILON = 1e-1;
/* ---------------------------------------------------------------------- */
ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr), rho0(nullptr),
Compute(lmp, narg, arg), chi(nullptr), fix_rheo(nullptr), rho0(nullptr),
norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr), compute_kernel(nullptr),
fix_pressure(nullptr)
@ -64,13 +62,12 @@ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
// between timesteps (fix property atom will handle callbacks)
int tmp1, tmp2;
int index = atom->find_custom("fp_store", tmp1, tmp2);
if (index == -1) {
index_fp_store = atom->find_custom("fp_store", tmp1, tmp2);
if (index_fp_store == -1) {
id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa));
index = atom->find_custom("fp_store", tmp1, tmp2);
index_fp_store = atom->find_custom("fp_store", tmp1, tmp2);
}
fp_store = atom->darray[index];
}
/* ---------------------------------------------------------------------- */
@ -114,7 +111,7 @@ void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr)
void ComputeRHEOInterface::compute_peratom()
{
int a, i, j, ii, jj, jnum, itype, fluidi, fluidj, status_match;
int a, i, j, ii, jj, jnum, fluidi, fluidj, status_match;
double xtmp, ytmp, ztmp, rsq, w, dot, dx[3];
int inum, *ilist, *jlist, *numneigh, **firstneigh;
@ -126,6 +123,7 @@ void ComputeRHEOInterface::compute_peratom()
int newton = force->newton;
int *status = atom->rheo_status;
double *rho = atom->rho;
double **fp_store = atom->darray[index_fp_store];
inum = list->inum;
ilist = list->ilist;
@ -151,7 +149,6 @@ void ComputeRHEOInterface::compute_peratom()
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
fluidi = !(status[i] & PHASECHECK);
jlist = firstneigh[i];
jnum = numneigh[i];
@ -212,9 +209,10 @@ void ComputeRHEOInterface::compute_peratom()
if (status[i] & PHASECHECK) {
if (normwf[i] != 0.0) {
// Stores rho for solid particles 1+Pw in Adami Adams 2012
rho[i] = MAX(EPSILON, fix_pressure->calc_rho(rho[i] / normwf[i], i));
// cap out at a tenth of equilibrium
rho[i] = MAX(0.1 * rho0[type[i]], fix_pressure->calc_rho(rho[i] / normwf[i], i));
} else {
rho[i] = rho0[itype];
rho[i] = rho0[type[i]];
}
}
}
@ -230,6 +228,7 @@ int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /
{
int m = 0;
double *rho = atom->rho;
double **fp_store = atom->darray[index_fp_store];
for (int i = 0; i < n; i++) {
int j = list[i];
@ -250,6 +249,8 @@ int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /
void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
{
double *rho = atom->rho;
double **fp_store = atom->darray[index_fp_store];
int m = 0;
int last = first + n;
for (int i = first; i < last; i++) {
@ -335,6 +336,7 @@ void ComputeRHEOInterface::store_forces()
double *mass = atom->mass;
double *rmass = atom->rmass;
double **f = atom->f;
double **fp_store = atom->darray[index_fp_store];
// When this is called, fp_store stores the pressure force
// After this method, fp_store instead stores non-pressure forces

3
src/RHEO/compute_rheo_interface.h
View File

@ -40,7 +40,8 @@ class ComputeRHEOInterface : public Compute {
double correct_rho(int);
void store_forces();
double *chi, **fp_store;
double *chi;
int index_fp_store;
class FixRHEO *fix_rheo;
private:

2
src/RHEO/compute_rheo_kernel.cpp
View File

@ -220,6 +220,8 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji);
if (kernel_style == QUINTIC)
return calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
if (kernel_style == RK0)
return calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
double wp;
int corrections_i = check_corrections(i);

2
src/RHEO/compute_rheo_property_atom.cpp
View File

@ -502,7 +502,7 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n)
void ComputeRHEOPropertyAtom::pack_nbond_shell(int n)
{
int *nbond = fix_oxidation->nbond;
int *nbond = atom->ivector[fix_oxidation->index_nb];
int *mask = atom->mask;
int nlocal = atom->nlocal;

19
src/RHEO/fix_rheo.cpp
View File

@ -39,6 +39,7 @@ using namespace RHEO_NS;
using namespace FixConst;
static const char cite_rheo[] =
"RHEO package: doi:10.1063/5.0228823\n\n"
"@article{Palermo2024,\n"
" journal = {Physics of Fluids},\n"
" title = {Reproducing hydrodynamics and elastic objects: A hybrid mesh-free model framework for dynamic multi-phase flows},\n"
@ -117,17 +118,17 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
shift_flag = 1;
memory->create(shift_type, n + 1, "rheo:shift_type");
for (i = 1; i <= n; i++) shift_type[i] = 1;
while (iarg < narg) { // optional sub-arguments
if (strcmp(arg[iarg], "scale/cross/type") == 0) {
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error);
while (iarg + 1 < narg) { // optional sub-arguments
if (strcmp(arg[iarg + 1], "scale/cross/type") == 0) {
if (iarg + 4 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error);
shift_cross_type_flag = 1;
shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
shift_cmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
shift_wmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
shift_scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
shift_wmin = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
iarg += 3;
} else if (strcmp(arg[iarg], "exclude/type") == 0) {
if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error);
utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error);
} else if (strcmp(arg[iarg + 1], "exclude/type") == 0) {
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error);
utils::bounds(FLERR, arg[iarg + 2], 1, n, nlo, nhi, error);
for (i = nlo; i <= nhi; i++) shift_type[i] = 0;
iarg += 1;
} else {

2
src/RHEO/fix_rheo_oxidation.cpp
View File

@ -38,6 +38,7 @@ using namespace FixConst;
enum { NONE, CONSTANT };
static const char cite_rheo_oxide[] =
"RHEO oxidation: doi:10.1016/j.apm.2024.02.027\n\n"
"@article{ApplMathModel.130.310,\n"
" title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
" journal = {Applied Mathematical Modelling},\n"
@ -109,7 +110,6 @@ void FixRHEOOxidation::init()
int tmp1, tmp2;
index_nb = atom->find_custom("shell_nbond", tmp1, tmp2);
if (index_nb == -1) error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation");
nbond = atom->ivector[index_nb];
// need a half neighbor list
auto req = neighbor->add_request(this, NeighConst::REQ_FULL);

4
src/RHEO/fix_rheo_oxidation.h
View File

@ -39,11 +39,11 @@ class FixRHEOOxidation : public Fix {
void post_force(int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int *nbond;
double rsurf, cut;
int index_nb;
private:
int btype, index_nb;
int btype;
double cutsq;
class NeighList *list;

5
src/RHEO/fix_rheo_thermal.cpp
View File

@ -221,11 +221,6 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
FixRHEOThermal::~FixRHEOThermal()
{
// Remove custom property if it exists
int tmp1, tmp2, index;
index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
if (index != -1) atom->remove_custom(index, 1, 0);
memory->destroy(cv_style);
memory->destroy(Tc_style);
memory->destroy(kappa_style);

6
src/RHEO/pair_rheo.cpp
View File

@ -115,7 +115,7 @@ void PairRHEO::compute(int eflag, int vflag)
double **fp_store, *chi;
if (compute_interface) {
fp_store = compute_interface->fp_store;
fp_store = atom->darray[compute_interface->index_fp_store];
chi = compute_interface->chi;
for (i = 0; i < atom->nmax; i++) {
@ -539,7 +539,7 @@ double PairRHEO::init_one(int i, int j)
int PairRHEO::pack_reverse_comm(int n, int first, double *buf)
{
double **fp_store = compute_interface->fp_store;
double **fp_store = atom->darray[compute_interface->index_fp_store];
int m = 0;
int last = first + n;
for (int i = first; i < last; i++) {
@ -555,7 +555,7 @@ int PairRHEO::pack_reverse_comm(int n, int first, double *buf)
void PairRHEO::unpack_reverse_comm(int n, int *list, double *buf)
{
double **fp_store = compute_interface->fp_store;
double **fp_store = atom->darray[compute_interface->index_fp_store];
int m = 0;
for (int i = 0; i < n; i++) {
int j = list[i];

4
src/comm_brick.cpp
View File

@ -209,8 +209,8 @@ void CommBrick::setup()
double cut = get_comm_cutoff();
if ((cut == 0.0) && (me == 0))
error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms "
"will be generated. Atoms may get lost.");
error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms will be generated. "
"Energies and forces may be wrong and atoms may get lost.");
if (mode == Comm::MULTI) {
double **cutcollectionsq = neighbor->cutcollectionsq;

4
src/comm_tiled.cpp
View File

@ -280,8 +280,8 @@ void CommTiled::setup()
double cut = get_comm_cutoff();
if ((cut == 0.0) && (me == 0))
error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms "
"will be generated. Atoms may get lost.");
error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms will be generated. "
"Energies and forces may be wrong and atoms may get lost.");
if (triclinic == 0) cutghost[0] = cutghost[1] = cutghost[2] = cut;
else {

2
src/displace_atoms.cpp
View File

@ -125,7 +125,7 @@ void DisplaceAtoms::command(int narg, char **arg)
else if (strcmp(arg[2],"y") == 0) d_dim = 1;
else if (strcmp(arg[2],"z") == 0) d_dim = 2;
else error->all(FLERR, 2, "Unknown displace_atoms ramp dimension {}", arg[2]);
if ((domain->dimension == 2) && (d_dim = 2))
if ((domain->dimension == 2) && (d_dim == 2))
error->all(FLERR, 2, "Must not displace atoms in z-direction with 2d system");
double d_lo,d_hi;

4
src/dump_custom.cpp
View File

@ -78,8 +78,8 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
int ioffset = 5;
expand = 0;
nfield = nargnew = utils::expand_args(FLERR,narg-5,&arg[5],1,earg,lmp);
if (earg != &arg[5]) expand = 1;
nfield = nargnew = utils::expand_args(FLERR,narg-ioffset,&arg[ioffset],1,earg,lmp);
if (earg != &arg[ioffset]) expand = 1;
// allocate field vectors

10
src/info.cpp
View File

@ -44,7 +44,9 @@
#include "region.h"
#include "update.h"
#include "variable.h"
#ifndef FMT_STATIC_THOUSANDS_SEPARATOR
#include "fmt/chrono.h"
#endif
#include <cctype>
#include <cmath>
@ -270,8 +272,16 @@ void Info::command(int narg, char **arg)
if (out == nullptr) return;
fputs("\nInfo-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n",out);
#if defined(FMT_STATIC_THOUSANDS_SEPARATOR)
{
time_t tv = time(nullptr);
struct tm *now = localtime(&tv);
utils::print(out, "Printed on {}", asctime(now));
}
#else
std::tm now = fmt::localtime(std::time(nullptr));
utils::print(out,"Printed on {}", std::asctime(&now));
#endif
if (flags & CONFIG) {
utils::print(out,"\nLAMMPS version: {} / {}\n", lmp->version, lmp->num_ver);

2
src/json.h
View File

@ -17,7 +17,7 @@
// wrapper around including the JSON parsing and writing class
// Do NOT include in any header file
#include "../third_party/nlohmann/json.hpp"
#include "nlohmann/json.hpp"
namespace LAMMPS_NS {
using json = ::nlohmann_lmp::json;

92
src/library.cpp
View File

@ -24,6 +24,7 @@
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "command.h"
#include "compute.h"
#include "domain.h"
#include "dump.h"
@ -42,6 +43,7 @@
#include "molecule.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "output.h"
#include "pair.h"
#if defined(LMP_PLUGIN)
@ -6119,6 +6121,96 @@ int lammps_find_compute_neighlist(void *handle, const char *id, int reqid) {
/* ---------------------------------------------------------------------- */
// helper Command class for a single neighbor list build
namespace LAMMPS_NS {
class NeighProxy : protected Command
{
public:
NeighProxy(class LAMMPS *lmp) : Command(lmp), neigh_idx(-1) {};
void command(int, char **) override;
int get_index() const { return neigh_idx; }
protected:
int neigh_idx;
};
}
void NeighProxy::command(int narg, char **arg)
{
neigh_idx = -1;
if (narg != 3) return;
auto *req = neighbor->add_request(this, arg[0]);
int flags = atoi(arg[1]);
double cutoff = atof(arg[2]);
req->apply_flags(flags);
if (cutoff > 0.0) req->set_cutoff(cutoff);
lmp->init();
// setup domain, communication and neighboring
// acquire ghosts and build standard neighbor lists
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
neighbor->build(1);
// build neighbor list this command needs based on earlier request
auto list = neighbor->find_list(this);
neighbor->build_one(list);
// find neigh list
for (int i = 0; i < neighbor->nlist; i++) {
NeighList *list = neighbor->lists[i];
if (this == list->requestor) {
neigh_idx = i;
break;
}
}
}
/** Build a single neighbor list in between runs and return its index
*
* A neighbor list request is created and the neighbor list built from a
* proxy command. The request ID is typically 0, but will be
* > 0 in case a compute has multiple neighbor list requests.
*
* \param handle pointer to a previously created LAMMPS instance cast to ``void *``.
* \param id Identifier of neighbor list request
* \param flags Neighbor list flags
* \param cutoff Custom neighbor list cutoff if > 0, use default cutoff if < 0
* \return return neighbor list index if valid, otherwise -1 */
int lammps_request_single_neighlist(void *handle, const char *id, int flags, double cutoff) {
auto lmp = (LAMMPS *)handle;
int idx = -1;
if (!lmp || !lmp->error || !lmp->neighbor) {
lammps_last_global_errormessage = fmt::format("ERROR: {}(): Invalid LAMMPS handle\n", FNERR);
return -1;
}
BEGIN_CAPTURE
{
NeighProxy proxy(lmp);
std::vector<std::string> args = {id, std::to_string(flags), std::to_string(cutoff)};
std::vector<char *> c_args;
std::transform(args.begin(), args.end(), std::back_inserter(c_args),
[](const std::string& s) { return (char *)s.c_str(); });
proxy.command(static_cast<int>(c_args.size()), c_args.data());
idx = proxy.get_index();
}
END_CAPTURE
return idx;
}
/* ---------------------------------------------------------------------- */
/** Return the number of entries in the neighbor list with given index
*
* \param handle pointer to a previously created LAMMPS instance cast to ``void *``.

12
src/library.h
View File

@ -109,6 +109,17 @@ enum _LMP_VAR_CONST {
LMP_VAR_STRING = 3 /*!< return value will be a string (catch-all) */
};
/** Neighbor list settings constants
*
* Must be kept in sync with the equivalent constants in ``python/lammps/constants.py``,
* ``fortran/lammps.f90``, ``tools/swig/lammps.i``, and
* ``examples/COUPLE/plugin/liblammpsplugin.h`` */
enum _LMP_NEIGH_CONST {
LMP_NEIGH_HALF = 0, /*!< request (default) half neighbor list */
LMP_NEIGH_FULL = 1, /*!< request full neighbor list */
};
/* Ifdefs to allow this file to be included in C and C++ programs */
#ifdef __cplusplus
@ -235,6 +246,7 @@ int lammps_create_atoms(void *handle, int n, const int64_t *id, const int *type,
int lammps_find_pair_neighlist(void *handle, const char *style, int exact, int nsub, int request);
int lammps_find_fix_neighlist(void *handle, const char *id, int request);
int lammps_find_compute_neighlist(void *handle, const char *id, int request);
int lammps_request_single_neighlist(void *handle, const char *id, int flags, double cutoff);
int lammps_neighlist_num_elements(void *handle, int idx);
void lammps_neighlist_element_neighbors(void *handle, int idx, int element, int *iatom,
int *numneigh, int **neighbors);

5
src/min.cpp
View File

@ -113,7 +113,10 @@ void Min::init()
// create fix needed for storing atom-based quantities
// will delete it at end of run
fix_minimize = dynamic_cast<FixMinimize *>(modify->add_fix("MINIMIZE all MINIMIZE"));
if (lmp->kokkos)
fix_minimize = dynamic_cast<FixMinimize *>(modify->add_fix("MINIMIZE all MINIMIZE/kk"));
else
fix_minimize = dynamic_cast<FixMinimize *>(modify->add_fix("MINIMIZE all MINIMIZE"));
// clear out extra global and per-atom dof
// will receive requests for new per-atom dof during pair init()

2831
third_party/nlohmann/json.hpp → src/nlohmann/json.hpp
View File

File diff suppressed because it is too large Load Diff

2
src/pair.cpp
View File

@ -90,6 +90,8 @@ Pair::Pair(LAMMPS *lmp) :
reinitflag = 1;
centroidstressflag = CENTROID_SAME;
atomic_energy_enable = 0;
// pair_modify settings
compute_flag = 1;

3
src/pair.h
View File

@ -54,6 +54,8 @@ class Pair : protected Pointers {
int single_enable; // 1 if single() routine exists
int born_matrix_enable; // 1 if born_matrix() routine exists
int single_hessian_enable; // 1 if single_hessian() routine exists
int atomic_energy_enable; // 1 if compute_atomic_energy() routine exists
int restartinfo; // 1 if pair style writes restart info
int respa_enable; // 1 if inner/middle/outer rRESPA routines
int one_coeff; // 1 if allows only one coeff * * call
@ -156,6 +158,7 @@ class Pair : protected Pointers {
virtual void compute_inner() {}
virtual void compute_middle() {}
virtual void compute_outer(int, int) {}
virtual double compute_atomic_energy(int, NeighList *) { return 0.0; }
virtual double single(int, int, int, int, double, double, double, double &fforce)
{

9
src/timer.cpp
View File

@ -15,7 +15,9 @@
#include "comm.h"
#include "error.h"
#ifndef FMT_STATIC_THOUSANDS_SEPARATOR
#include "fmt/chrono.h"
#endif
#include "tokenizer.h"
#include <cstring>
@ -264,8 +266,15 @@ void Timer::modify_params(int narg, char **arg)
// format timeout setting
std::string timeout = "off";
if (_timeout >= 0.0) {
#if defined(FMT_STATIC_THOUSANDS_SEPARATOR)
char outstr[200];
struct tm *tv = gmtime(&((time_t) _timeout));
strftime(outstr, 200, "%02d:%M:%S", tv);
timeout = outstr;
#else
std::tm tv = fmt::gmtime((std::time_t) _timeout);
timeout = fmt::format("{:02d}:{:%M:%S}", tv.tm_yday * 24 + tv.tm_hour, tv);
#endif
}
utils::logmesg(lmp, "New timer settings: style={} mode={} timeout={}\n", timer_style[_level],

9
src/utils.cpp
View File

@ -20,7 +20,9 @@
#include "domain.h"
#include "error.h"
#include "fix.h"
#ifndef FMT_STATIC_THOUSANDS_SEPARATOR
#include "fmt/chrono.h"
#endif
#include "input.h"
#include "label_map.h"
#include "memory.h"
@ -1953,8 +1955,15 @@ int utils::date2num(const std::string &date)
std::string utils::current_date()
{
time_t tv = time(nullptr);
#if defined(FMT_STATIC_THOUSANDS_SEPARATOR)
char outstr[200];
struct tm *today = localtime(&tv);
strftime(outstr, 200, "%Y-%m-%d", today);
return std::string(outstr);
#else
std::tm today = fmt::localtime(tv);
return fmt::format("{:%Y-%m-%d}", today);
#endif
}
/* ----------------------------------------------------------------------

14
third_party/README vendored
View File

@ -1 +1,13 @@
This folder contains copies of third-party software
This folder contains copies of third-party software.
These packages may be either included in the source
code tarball or downloaded on-the-fly during compilation
as needed. In some cases, some small customizations
are made to avoid namespace conflicts where the
included third party software needs workarounds for
issues on platforms (OS or compiler related) to
function consistently with LAMMPS.
The following software source packages are included
Folder: Package: Version: Original URL
nlohmann/json JSON for Modern C++ 3.12.0 https://github.com/nlohmann/json/releases/tag/v3.12.0

5
tools/coding_standard/README
View File

@ -6,3 +6,8 @@ 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
versiontags.py detects .. versionadded:: and .. versionchanged:: with pending version date
Example usage: python3 versiontags.py <DIRECTORY>
where <DIRECTORY> is the path of a LAMMPS directory containing subdirs src/, python/, doc/src/, etc.

17
tools/lammps-gui/flagwarnings.cpp
View File

@ -23,18 +23,22 @@
// workaround for Qt-5.12
#if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
namespace QColorConstants {
const QColor Red = QColor::fromRgb(0xff, 0x00, 0x00);
const QColor Red = QColor::fromRgb(0xff, 0x00, 0x00);
const QColor Blue = QColor::fromRgb(0x00, 0x00, 0xff);
} // namespace QColorConstants
#endif
FlagWarnings::FlagWarnings(QLabel *label, QTextDocument *parent) :
QSyntaxHighlighter(parent), isWarning(QStringLiteral("^(ERROR|WARNING).*$")), summary(label),
QSyntaxHighlighter(parent), isWarning(QStringLiteral("^(ERROR|WARNING).*$")),
isURL(QStringLiteral("^.*(https://docs.lammps.org/err[0-9]+).*$")), summary(label),
document(parent)
{
nwarnings = nlines = 0;
formatWarning.setForeground(QColorConstants::Red);
formatWarning.setFontWeight(QFont::Bold);
formatURL.setForeground(QColorConstants::Blue);
formatURL.setFontWeight(QFont::Bold);
}
void FlagWarnings::highlightBlock(const QString &text)
@ -42,11 +46,20 @@ void FlagWarnings::highlightBlock(const QString &text)
// nothing to do for empty lines
if (text.isEmpty()) return;
// highlight errors or warnings
auto match = isWarning.match(text);
if (match.hasMatch()) {
++nwarnings;
setFormat(match.capturedStart(0), match.capturedLength(0), formatWarning);
}
// highlight ErrorURL links
match = isURL.match(text);
if (match.hasMatch()) {
setFormat(match.capturedStart(1), match.capturedLength(1), formatURL);
}
// update error summary label
if (document && summary) {
summary->setText(
QString("%1 Warnings / Errors - %2 Lines").arg(nwarnings).arg(document->lineCount()));

2
tools/lammps-gui/flagwarnings.h
View File

@ -33,7 +33,9 @@ protected:
private:
QRegularExpression isWarning;
QRegularExpression isURL;
QTextCharFormat formatWarning;
QTextCharFormat formatURL;
QLabel *summary;
QTextDocument *document;
int nwarnings, nlines;

31
tools/lammps-gui/help_index.table
View File

@ -42,6 +42,7 @@ angle_lepton.html angle_style lepton
angle_lepton.html angle_style lepton/omp
angle_mesocnt.html angle_style mesocnt
angle_mm3.html angle_style mm3
angle_mwlc.html angle_style mwlc
angle_none.html angle_style none
angle_quartic.html angle_style quartic
angle_quartic.html angle_style quartic/omp
@ -129,6 +130,7 @@ atom_modify.html atom_modify
atom_style.html atom_style
balance.html balance
bond_bpm_rotational.html bond_style bpm/rotational
bond_bpm_spring_plastic.html bond_style bpm/spring/plastic
bond_bpm_spring.html bond_style bpm/spring
bond_class2.html bond_style class2
bond_class2.html bond_style class2/kk
@ -233,6 +235,8 @@ compute_event_displace.html compute event/displace
compute_fabric.html compute fabric
compute_fep.html compute fep
compute_fep_ta.html compute fep/ta
compute_gaussian_grid_local.html compute gaussian/grid/local
compute_gaussian_grid_local.html compute gaussian/grid/local/kk
compute_global_atom.html compute global/atom
compute_group_group.html compute group/group
compute_gyration_chunk.html compute gyration/chunk
@ -315,7 +319,9 @@ compute_smd_vol.html compute smd/vol
compute_sna_atom.html compute sna/atom
compute_sna_atom.html compute snad/atom
compute_sna_atom.html compute sna/grid
compute_sna_atom.html compute sna/grid/kk
compute_sna_atom.html compute sna/grid/local
compute_sna_atom.html compute sna/grid/local/kk
compute_sna_atom.html compute snap
compute_sna_atom.html compute snav/atom
compute_sph_e_atom.html compute sph/e/atom
@ -358,6 +364,7 @@ compute_temp_sphere.html compute temp/sphere
compute_temp_uef.html compute temp/uef
compute_ti.html compute ti
compute_torque_chunk.html compute torque/chunk
compute_vacf_chunk.html compute vacf/chunk
compute_vacf.html compute vacf
compute_vcm_chunk.html compute vcm/chunk
compute_viscosity_cos.html compute viscosity/cos
@ -396,6 +403,7 @@ dihedral_hybrid.html dihedral_style hybrid/kk
dihedral_lepton.html dihedral_style lepton
dihedral_lepton.html dihedral_style lepton/omp
dihedral_multi_harmonic.html dihedral_style multi/harmonic
dihedral_multi_harmonic.html dihedral_style multi/harmonic/kk
dihedral_multi_harmonic.html dihedral_style multi/harmonic/omp
dihedral_nharmonic.html dihedral_style nharmonic
dihedral_nharmonic.html dihedral_style nharmonic/omp
@ -436,6 +444,7 @@ dump.html dump custom
dump.html dump custom/gz
dump.html dump custom/zstd
dump.html dump dcd
dump.html dump extxyz
dump.html dump grid
dump.html dump grid/vtk
dump.html dump local
@ -487,7 +496,9 @@ fix_brownian.html fix brownian/asphere
fix_brownian.html fix brownian/sphere
fix_charge_regulation.html fix charge/regulation
fix_cmap.html fix cmap
fix_cmap.html fix cmap/kk
fix_colvars.html fix colvars
fix_colvars.html fix colvars/kk
fix_controller.html fix controller
fix_damping_cundall.html fix damping/cundall
fix_deform_pressure.html fix deform/pressure
@ -504,6 +515,7 @@ fix_drude_transform.html fix drude/transform/direct
fix_drude_transform.html fix drude/transform/inverse
fix_dt_reset.html fix dt/reset
fix_dt_reset.html fix dt/reset/kk
fix_efield_lepton.html fix efield/lepton
fix_efield.html fix efield
fix_efield.html fix efield/kk
fix_efield.html fix efield/tip4p
@ -614,6 +626,7 @@ fix_nve_dotc_langevin.html fix nve/dotc/langevin
fix_nve_dot.html fix nve/dot
fix_nve_eff.html fix nve/eff
fix_nve_limit.html fix nve/limit
fix_nve_limit.html fix nve/limit/kk
fix_nve_line.html fix nve/line
fix_nve_manifold_rattle.html fix nve/manifold/rattle
fix_nve_noforce.html fix nve/noforce
@ -647,7 +660,9 @@ fix_pafi.html fix pafi
fix_pair.html fix pair
fix_phonon.html fix phonon
fix_pimd.html fix pimd/langevin
fix_pimd.html fix pimd/langevin/bosonic
fix_pimd.html fix pimd/nvt
fix_pimd.html fix pimd/nvt/bosonic
fix_planeforce.html fix planeforce
fix_plumed.html fix plumed
fix_poems.html fix poems
@ -670,6 +685,8 @@ fix_qeq_comb.html fix qeq/comb/omp
fix_qeq_reaxff.html fix qeq/reaxff
fix_qeq_reaxff.html fix qeq/reaxff/kk
fix_qeq_reaxff.html fix qeq/reaxff/omp
fix_qeq_rel_reaxff.html fix qeq/rel/reaxff
fix_qeq.html fix qeq/ctip
fix_qeq.html fix qeq/dynamic
fix_qeq.html fix qeq/fire
fix_qeq.html fix qeq/point
@ -677,11 +694,13 @@ fix_qeq.html fix qeq/shielded
fix_qeq.html fix qeq/slater
fix_qmmm.html fix qmmm
fix_qtb.html fix qtb
fix_qtpie_reaxff.html fix qtpie/reaxff
fix_reaxff_bonds.html fix reaxff/bonds
fix_reaxff_bonds.html fix reaxff/bonds/kk
fix_reaxff_species.html fix reaxff/species
fix_reaxff_species.html fix reaxff/species/kk
fix_recenter.html fix recenter
fix_recenter.html fix recenter/kk
fix_restrain.html fix restrain
fix_rheo_oxidation.html fix rheo/oxidation
fix_rheo_pressure.html fix rheo/pressure
@ -772,6 +791,7 @@ fix_wall_reflect.html fix wall/reflect
fix_wall_reflect.html fix wall/reflect/kk
fix_wall_reflect_stochastic.html fix wall/reflect/stochastic
fix_wall_region.html fix wall/region
fix_wall_region.html fix wall/region/kk
fix_wall.html fix wall/colloid
fix_wall.html fix wall/harmonic
fix_wall.html fix wall/lepton
@ -783,6 +803,7 @@ fix_wall.html fix wall/morse
fix_wall.html fix wall/table
fix_wall_srd.html fix wall/srd
fix_widom.html fix widom
geturl.html geturl
group2ndx.html group2ndx
group2ndx.html ndx2group
group.html group
@ -919,6 +940,7 @@ pair_born.html pair_style born/gpu
pair_born.html pair_style born/omp
pair_bpm_spring.html pair_style bpm/spring
pair_brownian.html pair_style brownian
pair_brownian.html pair_style brownian/kk
pair_brownian.html pair_style brownian/omp
pair_brownian.html pair_style brownian/poly
pair_brownian.html pair_style brownian/poly/omp
@ -983,6 +1005,7 @@ pair_comb.html pair_style comb/omp
pair_cosine_squared.html pair_style cosine/squared
pair_coul_diel.html pair_style coul/diel
pair_coul_diel.html pair_style coul/diel/omp
pair_coul.html pair_style coul/ctip
pair_coul.html pair_style coul/cut
pair_coul.html pair_style coul/cut/global
pair_coul.html pair_style coul/cut/global/omp
@ -1049,6 +1072,7 @@ pair_dipole.html pair_style lj/long/dipole/long
pair_dipole.html pair_style lj/sf/dipole/sf
pair_dipole.html pair_style lj/sf/dipole/sf/gpu
pair_dipole.html pair_style lj/sf/dipole/sf/omp
pair_dispersion_d3.html pair_style dispersion/d3
pair_dpd_coul_slater_long.html pair_style dpd/coul/slater/long
pair_dpd_coul_slater_long.html pair_style dpd/coul/slater/long/gpu
pair_dpd_ext.html pair_style dpd/ext
@ -1153,8 +1177,12 @@ pair_gw.html pair_style gw/zbl
pair_harmonic_cut.html pair_style harmonic/cut
pair_harmonic_cut.html pair_style harmonic/cut/omp
pair_hbond_dreiding.html pair_style hbond/dreiding/lj
pair_hbond_dreiding.html pair_style hbond/dreiding/lj/angleoffset
pair_hbond_dreiding.html pair_style hbond/dreiding/lj/angleoffset/omp
pair_hbond_dreiding.html pair_style hbond/dreiding/lj/omp
pair_hbond_dreiding.html pair_style hbond/dreiding/morse
pair_hbond_dreiding.html pair_style hbond/dreiding/morse/angleoffset
pair_hbond_dreiding.html pair_style hbond/dreiding/morse/angleoffset/omp
pair_hbond_dreiding.html pair_style hbond/dreiding/morse/omp
pair_hdnnp.html pair_style hdnnp
pair_hybrid.html pair_style hybrid
@ -1236,6 +1264,8 @@ pair_lj_long.html pair_style lj/long/coul/long/omp
pair_lj_long.html pair_style lj/long/coul/long/opt
pair_lj_long.html pair_style lj/long/tip4p/long
pair_lj_long.html pair_style lj/long/tip4p/long/omp
pair_lj_pirani.html pair_style lj/pirani
pair_lj_pirani.html pair_style lj/pirani/omp
pair_lj_relres.html pair_style lj/relres
pair_lj_relres.html pair_style lj/relres/omp
pair_lj.html pair_style lj/cut
@ -1469,6 +1499,7 @@ quit.html quit
read_data.html read_data
read_dump.html read_dump
read_restart.html read_restart
region2vmd.html region2vmd
region.html region
replicate.html replicate
rerun.html rerun

3
tools/lammps-gui/lammps-gui.appdata.xml
View File

@ -62,7 +62,10 @@
Add option to visualize molecules defined through the molecule command
Add text fields for editing plot title and axis labels for charts
Add option to automatically open tutorial websites (enabled by default)
Add citation reminder comment to empty documents and the About dialog currently pointing to arXiv.
Add preferences tab for charts to set default for title, plot colors, smooth/raw plot, smooth params
Highlight error URLs with pointers to additional explanations in log window
Double-click on highlighted URL opens it in web browser. Also available via context menu.
</description>
</release>
<release version="1.6.12" timestamp="1734890080">

9
tools/lammps-gui/lammpsgui.cpp
View File

@ -68,6 +68,8 @@
static const QString blank(" ");
static constexpr int BUFLEN = 256;
static const QString citeme("# When using LAMMPS-GUI in your project, please cite: "
"https://arxiv.org/abs/2503.14020\n");
LammpsGui::LammpsGui(QWidget *parent, const QString &filename) :
QMainWindow(parent), ui(new Ui::LammpsGui), highlighter(nullptr), capturer(nullptr),
@ -78,6 +80,8 @@ LammpsGui::LammpsGui(QWidget *parent, const QString &filename) :
{
docver = "";
ui->setupUi(this);
ui->textEdit->document()->setPlainText(citeme);
ui->textEdit->document()->setModified(false);
this->setCentralWidget(ui->textEdit);
highlighter = new Highlighter(ui->textEdit->document());
capturer = new StdCapture;
@ -413,7 +417,8 @@ LammpsGui::~LammpsGui()
void LammpsGui::new_document()
{
current_file.clear();
ui->textEdit->document()->setPlainText(QString());
ui->textEdit->document()->setPlainText(citeme);
ui->textEdit->document()->setModified(false);
if (lammps.is_running()) {
stop_run();
@ -1487,10 +1492,10 @@ void LammpsGui::about()
info = std::string(info, start, end - start);
}
info += citeme.toStdString();
to_clipboard += info.c_str();
#if QT_CONFIG(clipboard)
QGuiApplication::clipboard()->setText(to_clipboard);
info += "(Note: this text has been copied to the clipboard)\n";
#endif
QMessageBox msg;

56
tools/lammps-gui/logwindow.cpp
View File

@ -18,6 +18,7 @@
#include <QAction>
#include <QApplication>
#include <QDesktopServices>
#include <QFile>
#include <QFileDialog>
#include <QGridLayout>
@ -37,6 +38,7 @@
const QString LogWindow::yaml_regex =
QStringLiteral("^(keywords:.*$|data:$|---$|\\.\\.\\.$| - \\[.*\\]$)");
const QString LogWindow::url_regex = QStringLiteral("^.*(https://docs.lammps.org/err[0-9]+).*$");
LogWindow::LogWindow(const QString &_filename, QWidget *parent) :
QPlainTextEdit(parent), filename(_filename), warnings(nullptr)
@ -198,8 +200,53 @@ void LogWindow::extract_yaml()
file.close();
}
void LogWindow::open_errorurl()
{
if (!errorurl.isEmpty()) QDesktopServices::openUrl(QUrl(errorurl));
}
void LogWindow::mouseDoubleClickEvent(QMouseEvent *event)
{
if (event->button() == Qt::LeftButton) {
// select the entire word (non-space text) under the cursor
// we need to do it in this complicated way, since QTextCursor does not recognize
// special characters as part of a word.
auto cursor = textCursor();
auto line = cursor.block().text();
int begin = qMin(cursor.positionInBlock(), line.length() - 1);
while (begin >= 0) {
if (line[begin].isSpace()) break;
--begin;
}
int end = begin + 1;
while (end < line.length()) {
if (line[end].isSpace()) break;
++end;
}
cursor.setPosition(cursor.position() - cursor.positionInBlock() + begin + 1);
cursor.movePosition(QTextCursor::NextCharacter, QTextCursor::KeepAnchor, end - begin - 1);
auto text = cursor.selectedText();
auto url = QRegularExpression(url_regex).match(text);
if (url.hasMatch()) {
errorurl = url.captured(1);
if (!errorurl.isEmpty()) {
QDesktopServices::openUrl(QUrl(errorurl));
return;
}
}
}
// forward event to parent class for all unhandled cases
QPlainTextEdit::mouseDoubleClickEvent(event);
}
void LogWindow::contextMenuEvent(QContextMenuEvent *event)
{
// reposition the cursor here, but only if there is no active selection
if (!textCursor().hasSelection()) setTextCursor(cursorForPosition(event->pos()));
// show augmented context menu
auto *menu = createStandardContextMenu();
menu->addSeparator();
@ -214,6 +261,15 @@ void LogWindow::contextMenuEvent(QContextMenuEvent *event)
action->setShortcut(QKeySequence(Qt::CTRL | Qt::Key_Y));
connect(action, &QAction::triggered, this, &LogWindow::extract_yaml);
}
// process line of text where the cursor is
auto text = textCursor().block().text().replace('\t', ' ').trimmed();
auto url = QRegularExpression(url_regex).match(text);
if (url.hasMatch()) {
errorurl = url.captured(1);
action = menu->addAction("Open &URL in Web Browser", this, &LogWindow::open_errorurl);
action->setIcon(QIcon(":/icons/help-browser.png"));
}
action = menu->addAction("&Jump to next warning or error", this, &LogWindow::next_warning);
action->setIcon(QIcon(":/icons/warning.png"));
action->setShortcut(QKeySequence(Qt::CTRL | Qt::Key_N));

4
tools/lammps-gui/logwindow.h
View File

@ -32,16 +32,20 @@ private slots:
void save_as();
void stop_run();
void next_warning();
void open_errorurl();
protected:
void closeEvent(QCloseEvent *event) override;
void mouseDoubleClickEvent(QMouseEvent *event) override;
void contextMenuEvent(QContextMenuEvent *event) override;
bool eventFilter(QObject *watched, QEvent *event) override;
bool check_yaml();
private:
QString filename;
QString errorurl;
static const QString yaml_regex;
static const QString url_regex;
FlagWarnings *warnings;
QLabel *summary;
};

2
tools/lammps-gui/org.lammps.lammps-gui.yml
View File

@ -1,6 +1,6 @@
id: org.lammps.lammps-gui
runtime: org.kde.Platform
runtime-version: "5.15-23.08"
runtime-version: "5.15-24.08"
sdk: org.kde.Sdk
command: lammps-gui
finish-args:

76
tools/swig/lammps.i
View File

@ -95,6 +95,17 @@ enum _LMP_VAR_CONST {
LMP_VAR_STRING = 3 /*!< return value will be a string (catch-all) */
};
/** Neighbor list settings constants
*
* Must be kept in sync with the equivalent constants in ``python/lammps/constants.py``,
* ``fortran/lammps.f90``, ``tools/swig/lammps.i``, and
* ``examples/COUPLE/plugin/liblammpsplugin.h`` */
enum _LMP_NEIGH_CONST {
LMP_NEIGH_HALF = 0, /*!< request (default) half neighbor list */
LMP_NEIGH_FULL = 1, /*!< request full neighbor list */
};
/*
extern void *lammps_open(int argc, char **argv, MPI_Comm comm, void **ptr);
*/
@ -143,32 +154,34 @@ extern int lammps_set_string_variable(void *, const char *, const char *);
extern int lammps_set_internal_variable(void *, const char *, double);
extern int lammps_variable_info(void *handle, int idx, char *buf, int bufsize);
extern double lammps_eval(void *handle, const char *expr);
extern void lammps_clearstep_compute(void *handle);
extern void lammps_addstep_compute(void *handle, void *nstep);
extern void lammps_addstep_compute_all(void *handle, void *nstep);
extern void lammps_gather_atoms(void *, char *, int, int, void *);
extern void lammps_gather_atoms_concat(void *, char *, int, int, void *);
extern void lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *);
extern void lammps_scatter_atoms(void *, char *, int, int, void *);
extern void lammps_scatter_atoms_subset(void *, char *, int, int, int, int *, void *);
extern void lammps_gather_atoms(void *, const char *, int, int, void *);
extern void lammps_gather_atoms_concat(void *, const char *, int, int, void *);
extern void lammps_gather_atoms_subset(void *, const char *, int, int, int, int *, void *);
extern void lammps_scatter_atoms(void *, const char *, int, int, void *);
extern void lammps_scatter_atoms_subset(void *, const char *, int, int, int, int *, void *);
extern void lammps_gather_bonds(void *handle, void *data);
extern void lammps_gather_angles(void *handle, void *data);
extern void lammps_gather_dihedrals(void *handle, void *data);
extern void lammps_gather_impropers(void *handle, void *data);
extern void lammps_gather(void *, char *, int, int, void *);
extern void lammps_gather_concat(void *, char *, int, int, void *);
extern void lammps_gather_subset(void *, char *, int, int, int, int *, void *);
extern void lammps_scatter(void *, char *, int, int, void *);
extern void lammps_scatter_subset(void *, char *, int, int, int, int *, void *);
extern int lammps_create_atoms(void *handle, int n, int *id, int *type,
double *x, double *v, int *image, int bexpand);
extern void lammps_gather(void *, const char *, int, int, void *);
extern void lammps_gather_concat(void *, const char *, int, int, void *);
extern void lammps_gather_subset(void *, const char *, int, int, int, int *, void *);
extern void lammps_scatter(void *, const char *, int, int, void *);
extern void lammps_scatter_subset(void *, const char *, int, int, int, int *, void *);
extern int lammps_create_atoms(void *handle, int n, const int *id, const int *type,
const double *x, const double *v, const int *image, int bexpand);
/*
extern int lammps_create_atoms(void *handle, int n, int64_t *id, int *type, */
extern int lammps_find_pair_neighlist(void*, char *, int, int, int);
extern int lammps_find_fix_neighlist(void*, char *, int);
extern int lammps_find_compute_neighlist(void*, char *, int);
extern int lammps_neighlist_num_elements(void*, int);
extern int lammps_find_pair_neighlist(void *, const char *, int, int, int);
extern int lammps_find_fix_neighlist(void *, const char *, int);
extern int lammps_find_compute_neighlist(void *, const char *, int);
extern int lammps_request_single_neighlist(void *, const char *, int, double);
extern int lammps_neighlist_num_elements(void *, int);
extern void lammps_neighlist_element_neighbors(void *, int, int, int *, int *, int ** );
extern int lammps_version(void *handle);
@ -344,27 +357,28 @@ extern void lammps_clearstep_compute(void *handle);
extern void lammps_addstep_compute(void *handle, void *nstep);
extern void lammps_addstep_compute_all(void *handle, void *nstep);
extern void lammps_gather_atoms(void *, char *, int, int, void *);
extern void lammps_gather_atoms_concat(void *, char *, int, int, void *);
extern void lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *);
extern void lammps_scatter_atoms(void *, char *, int, int, void *);
extern void lammps_scatter_atoms_subset(void *, char *, int, int, int, int *, void *);
extern void lammps_gather_atoms(void *, const char *, int, int, void *);
extern void lammps_gather_atoms_concat(void *, const char *, int, int, void *);
extern void lammps_gather_atoms_subset(void *, const char *, int, int, int, int *, void *);
extern void lammps_scatter_atoms(void *, const char *, int, int, void *);
extern void lammps_scatter_atoms_subset(void *, const char *, int, int, int, int *, void *);
extern void lammps_gather_bonds(void *handle, void *data);
extern void lammps_gather_angles(void *handle, void *data);
extern void lammps_gather_dihedrals(void *handle, void *data);
extern void lammps_gather_impropers(void *handle, void *data);
extern void lammps_gather(void *, char *, int, int, void *);
extern void lammps_gather_concat(void *, char *, int, int, void *);
extern void lammps_gather_subset(void *, char *, int, int, int, int *, void *);
extern void lammps_scatter(void *, char *, int, int, void *);
extern void lammps_scatter_subset(void *, char *, int, int, int, int *, void *);
extern int lammps_create_atoms(void *handle, int n, int *id, int *type,
double *x, double *v, int *image, int bexpand);
extern void lammps_gather(void *, const char *, int, int, void *);
extern void lammps_gather_concat(void *, const char *, int, int, void *);
extern void lammps_gather_subset(void *, const char *, int, int, int, int *, void *);
extern void lammps_scatter(void *, const char *, int, int, void *);
extern void lammps_scatter_subset(void *, const char *, int, int, int, int *, void *);
extern int lammps_create_atoms(void *handle, int n, const int *id, const int *type,
const double *x, const double *v, const int *image, int bexpand);
/*
extern int lammps_create_atoms(void *handle, int n, int64_t *id, int *type, */
extern int lammps_find_pair_neighlist(void*, char *, int, int, int);
extern int lammps_find_fix_neighlist(void*, char *, int);
extern int lammps_find_compute_neighlist(void*, char *, int);
extern int lammps_find_pair_neighlist(void*, const char *, int, int, int);
extern int lammps_find_fix_neighlist(void*, const char *, int);
extern int lammps_find_compute_neighlist(void*, const char *, int);
extern int lammps_request_single_neighlist(void *, const char *, int, double);
extern int lammps_neighlist_num_elements(void*, int);
extern void lammps_neighlist_element_neighbors(void *, int, int, int *, int *, int ** );

1
tools/tabulate/README.md
View File

@ -27,6 +27,7 @@ Please see the individual tabulation scripts in this folder for examples:
| wall_harmonic_tabulate.py | creates a table for fix wall/table with a simple repulsive harmonic potential |
| wall_multi_tabulate.py | creates a table for fix wall/table with multiple tables |
| pair_bi_tabulate.py | creates a table from radial distribution file using Boltzmann Inversion |
| plot_forces.py | plots and extracts tabulated forces from table files |
Common command line flags:

12
tools/tabulate/angle_harmonic_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import AngleTabulate

12
tools/tabulate/bond_morse_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import BondTabulate

12
tools/tabulate/dihedral_harmonic_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import DihedralTabulate

12
tools/tabulate/pair_bi_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python3
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import PairTabulate
import sys

12
tools/tabulate/pair_hybrid_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import PairTabulate

12
tools/tabulate/pair_lj_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import PairTabulate

12
tools/tabulate/pair_zbladd_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import PairTabulate

290
tools/tabulate/plot_forces.py Executable file
View File

@ -0,0 +1,290 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
# Author: Germain Clavier (Unicaen), germain.clavier at unicaen.fr
"""
Plot LAMMPS tabulated forces.
"""
import argparse
import numpy as np
import os
import logging
import sys
from matplotlib import pyplot as plt
logger = logging.getLogger(__name__)
units = {
"lj": {
"Distance": "Reduced units",
"Energy": "Reduced units",
"Force": "Reduced units",
"kb": 1,
},
"real": {
"Distance": "[A]",
"Energy": "[kcal/mol]",
"Force": "[kcal/mol/A]",
"kb": 0.001985875,
},
"metal": {
"Distance": "[A]",
"Energy": "[eV]",
"Force": "[eV/A]",
"kb": 8.6173332e-5,
},
"si": {"Distance": "[m]", "Energy": "[J]", "Force": "[N]", "kb": 1.380649e-23},
"cgs": {
"Distance": "[cm]",
"Energy": "[ergs]",
"Force": "[dynes]",
"kb": 1.3806504e-16,
},
"electron": {
"Distance": "[Bohr]",
"Energy": "[Hartrees]",
"Force": "[Hartree/Bohr]",
"kb": 3.16681534e-6,
},
"micro": {
"Distance": "[um]",
"Energy": "[pg·um^2/us^2]",
"Force": "[pg·um/us^2]",
"kb": 1.3806504e-8,
},
"nano": {
"Distance": "[nm]",
"Energy": "[ag·nm^2/ns^2]",
"Force": "[ag·nm/ns^2]",
"kb": 0.013806504,
},
}
def compute_energy(tp):
r = tp[0]
fo = tp[2]
e = np.zeros(r.shape)
for i, (ri, fi) in enumerate(zip(r, fo)):
if i == 0:
continue
dr = ri - r[i - 1]
e[i] = e[i - 1] - dr * fo[i - 1]
e -= e[-1]
return e
def main():
parser = argparse.ArgumentParser(
description="""
Plots LAMMPS tabulated forces. This script takes a table
file as an input and plots all the tabulated forces inside with their
corresponding energy. The forces label is the token used to name the
force in the file. It can be used to output all the forces in separate
files and/or recompute the energy from forces through finite difference
(assuming e(rc)=0). This script requires the matplotlib and numpy
Python libraries. Bitmap format is not supported.
"""
)
parser.add_argument(
"-u",
"--units",
dest="units",
default="real",
help="Units of the file (LAMMPS units system)",
)
parser.add_argument(
"-f",
"--file",
dest="infile",
default="",
help="File to read",
)
parser.add_argument(
"-x",
dest="xrange",
default="",
help="xrange separated by : (for negative values use the '=' sign: -x=-3:10)",
)
parser.add_argument(
"-y",
dest="yrange",
default="",
help="yrange separated by :",
)
parser.add_argument(
"-t",
dest="temp",
default=None,
type=float,
help="temperature for KbT plot [default none]",
)
parser.add_argument(
"-d",
"--diff-num",
dest="recompute",
action="store_true",
help="Recompute the energies from forces and distances through finite differences",
)
parser.add_argument(
"-e",
dest="extract",
action="store_true",
help="Extract the forces in separate files",
)
args = parser.parse_args()
logging.basicConfig(level=logging.INFO)
##########
# Manage arguments
udistance = units[args.units]["Distance"]
uenergy = units[args.units]["Energy"]
uforce = units[args.units]["Force"]
kb = units[args.units]["kb"]
rlabel = " ".join(["Rij", udistance])
elabel = " ".join(["E", uenergy])
flabel = " ".join(["F", uforce])
etitle = "Energy"
ftitle = "Force"
font = "DejaVu Sans"
fontsize = 30
infile = args.infile
if not os.path.isfile(infile):
logger.error("Input file not found")
sys.exit(1)
toplot = []
with open(infile, "r") as f:
lines = iter(f.readlines())
while True:
try:
r = []
force = []
ener = []
tok = []
while not tok:
tok = next(lines).partition("#")[0].rstrip()
logger.info("Found {} token".format(tok))
infos = next(lines).split()
npoints = int(infos[1])
next(lines)
if "bitmap" in infos:
logger.info("Unsupported bitmap format for token {:s}".format(tok))
for _ in range(npoints):
continue
else:
for i in range(npoints):
line = next(lines).split()
r.append(float(line[1]))
ener.append(float(line[2]))
force.append(float(line[3]))
r = np.array(r)
ener = np.array(ener)
force = np.array(force)
toplot.append([r, ener, force, tok])
tok = []
next(lines)
except StopIteration:
break
if args.recompute:
etitle = "Estimated energy"
for tp in toplot:
tp[1] = compute_energy(tp)
fig, axes = plt.subplots(1, 2)
for tp in toplot:
axes[0].plot(tp[0], tp[1], label=tp[3], linewidth=3)
axes[1].plot(tp[0], tp[2], label=tp[3], linewidth=3)
hmin, hmax = axes[1].get_xlim()
axes[1].hlines(0, hmin, hmax, color="black", linewidth=3, linestyles="dashdot")
if args.temp:
if args.temp > 0:
hmin, hmax = axes[0].get_xlim()
axes[0].hlines(
kb * args.temp,
hmin,
hmax,
color="orange",
label=r"$k_BT$",
linewidth=3,
linestyles="dashdot",
)
axes[0].text(hmax / 2.0, kb * args.temp, "KbT", fontsize=0.7 * fontsize)
logger.info("KbT value= {:e} {:s}".format(kb * args.temp, uenergy))
else:
logger.info("Invalid temperature value: {:e}".format(args.temp))
if args.xrange:
xmin, xmax = list(map(float, args.xrange.split(":")))
axes[0].set_xlim(xmin, xmax)
axes[1].set_xlim(xmin, xmax)
if args.yrange:
ymin, ymax = list(map(float, args.yrange.split(":")))
axes[0].set_ylim(ymin, ymax)
axes[1].set_ylim(ymin, ymax)
# Setting axes 0
axes[0].set_title(etitle, fontsize=fontsize)
axes[0].set_xlabel(
rlabel, fontname=font, fontsize=fontsize
) # xlabel name, size 30pts
axes[0].set_ylabel(
elabel, fontname=font, fontsize=fontsize
) # ylabel name, size 30pts
axes[0].tick_params(
axis="both", which="major", labelsize=fontsize
) # Biggers ticks, bigger tick labels!
# Setting axes 1
axes[1].set_title(ftitle, fontsize=fontsize)
axes[1].legend(frameon=False, fontsize=fontsize) # Fat font, no frame
axes[1].set_xlabel(
rlabel, fontname=font, fontsize=fontsize
)
axes[1].set_ylabel(
flabel, fontname=font, fontsize=fontsize
)
axes[1].tick_params(
axis="both", which="major", labelsize=fontsize
)
figManager = plt.get_current_fig_manager()
figManager.window.showMaximized()
plt.show()
if args.extract:
for tp in toplot:
outfile = "".join([tp[3], ".plot"])
logger.info("Writing file {}".format(outfile))
with open(outfile, "w") as f:
f.write("# {} force extracted from {}\n".format(tp[3], infile))
f.write("# {:^20} {:^20} {:^20}\n".format("r", "energy", "force"))
for a, b, c in zip(tp[0], tp[1], tp[2]):
f.write("{:>18.16e} {:>18.16e} {:>18.16e}\n".format(a, b, c))
return
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
raise SystemExit("User interruption.")

12
tools/tabulate/wall_harmonic_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import WallTabulate

12
tools/tabulate/wall_multi_tabulate.py
View File

@ -1,4 +1,16 @@
#!/usr/bin/env python
# ----------------------------------------------------------------------
# LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
# https://www.lammps.org/ Sandia National Laboratories
# LAMMPS Development team: developers@lammps.org
#
# Copyright (2003) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
#
# See the README file in the top-level LAMMPS directory.
# -------------------------------------------------------------------------
from tabulate import WallTabulate
import os, sys

2
unittest/commands/test_kim_commands.cpp
View File

@ -636,7 +636,7 @@ TEST_F(KimCommandsTest, kim_query)
"temperature_units=[K]");
END_HIDE_OUTPUT();
ASSERT_THAT(variable->retrieve("alpha"), StrEq("1.654960564704273e-05"));
ASSERT_THAT(variable->retrieve("alpha"), StrEq("1.656579473023212e-05"));
BEGIN_HIDE_OUTPUT();
command("clear");

274
unittest/force-styles/tests/atomic-pair-buck_coul_cut_qeq_point.yaml
View File

@ -1,7 +1,7 @@
---
lammps_version: 17 Feb 2022
lammps_version: 2 Apr 2025
tags: unstable
date_generated: Fri Mar 18 22:17:36 2022
date_generated: Fri Apr 25 19:44:59 2025
epsilon: 7.5e-13
skip_tests: gpu intel single
prerequisites: ! |
@ -21,7 +21,7 @@ pre_commands: ! |
displace_atoms all random 0.1 0.1 0.1 623426
mass 1 28.0
mass 2 16.0
set type 1 type/fraction 2 0.666667 998877
set type 1 type/ratio 2 0.666667 998877
set type 1 charge 0.8
set type 2 charge 0.4
velocity all create 100 4534624 loop geom
@ -35,142 +35,142 @@ pair_coeff: ! |
2 2 1388.77 0.3623188 175.0
extract: ! ""
natoms: 64
init_vdwl: 106.41860675757
init_coul: -129.9140700924124
init_vdwl: -24.42807736156882
init_coul: -131.75120059001813
init_stress: ! |-
-2.8069812877025680e+02 -4.4806768911794171e+02 -5.0378472147221044e+02 3.9357896002405380e+02 -4.6782123097483543e+02 1.9690042286173212e+02
-6.8611278946113976e+02 -8.4852099966145420e+02 -7.2493943802866158e+02 8.7199041253461814e+02 -6.9309149754481655e+02 1.0108767842057036e+03
init_forces: ! |2
1 7.9125122896296922e+00 -2.1014839690852849e+00 1.0994080971957654e+01
2 -8.6188789139298905e+00 -1.8333074427098923e+00 6.9043233196778493e+00
3 7.8498612200411868e+00 -2.1754399089581845e+01 -2.8862198492435933e+00
4 6.0059474900265855e+00 2.0200604615492711e+01 -3.6502720655785064e+00
5 -1.0414088486908696e+01 -1.2032673161543727e+01 -8.1683250459138641e+00
6 2.4526971509997360e+00 1.6151940503875748e+00 3.9699353225387135e+00
7 1.1484521325651471e+01 2.6251346011917382e+00 -5.0509149729671243e-01
8 -8.6634039445313658e+00 -8.3084485039834881e+00 3.9070922142448175e+00
9 7.9495750681929742e+00 -9.0507972374486361e+00 -4.4437699190231843e+00
10 -1.1184211738488107e+01 -2.0985043628368132e+00 8.0230797778549956e+00
11 2.8166289803609963e+01 -2.6396601038400636e+01 -7.0398077372607645e+01
12 1.8384817166363799e+01 4.9304421068745743e+00 -1.8773815008331542e+01
13 -8.5901832560112172e+00 -1.1830174021232128e+01 1.0687567118271307e+01
14 -2.7556923264654802e+00 1.7745847764979608e+00 -9.2066420206368988e+00
15 4.6735302548764643e+00 -3.2709618690940427e+00 -2.0693142194244816e+00
16 -6.0547965516444906e+00 -2.2600378957272116e+00 3.2677368381468008e-01
17 -9.1378380656225406e+00 1.1908878845184031e+01 -1.7535407559841292e+00
18 -1.2710024935571287e+00 -1.8390249623807247e+00 -5.0661022351703764e+00
19 4.7508234682845059e+00 -4.7660218202653398e-01 1.9578845982916517e+00
20 -2.4044494501312546e+00 1.3188349248766018e+01 -4.2723195450669484e+00
21 2.5545939677742822e+00 -7.6599926401155196e+00 -1.9913355105865760e+00
22 -1.2721065348405212e+01 2.2605336717084196e+01 -4.7468563523008402e+00
23 2.4611421112071081e+00 -1.4649327715865049e-01 3.5863258472083799e+00
24 -4.3893602786163104e+00 2.1969519282498890e+01 4.4720052177174487e+00
25 2.0238917888202160e+01 -1.5707849435118018e+01 5.7119813878581205e+00
26 -1.2244546012542273e+01 8.7637388770168165e+00 1.1554289778536326e+01
27 2.6315583772500917e+00 5.2366381109035203e+00 5.6567854473912291e+00
28 1.6069342319200945e+01 2.2379145739762505e+00 1.2146106338402367e+01
29 1.6156898389873891e+01 -1.3554238523505591e+01 -9.1896135574796052e+00
30 -8.4056006337974853e+00 7.0469844075145707e+00 6.3110771274676658e+00
31 -1.3273848195602289e-01 -6.8307139969335129e+00 7.6340671688916251e-01
32 -3.5499992242635152e+01 1.6857356356036899e+01 5.1954595281751772e+01
33 -1.5976858206617210e+01 -2.1986638912202295e+01 1.5916728791281200e+01
34 9.9892552480335084e+00 2.4910832782965171e+00 -3.5317186985424054e+00
35 7.0710369156676450e+01 -3.5579041343629925e+01 -5.3907803084513858e+01
36 -5.7367330604749455e+01 -7.0169041999492791e+01 7.3977418913261829e+01
37 5.5359814356241074e+01 7.2104773934571938e+01 -6.4062698737301162e+01
38 7.3137814516005992e+00 -6.9993845661823428e+00 -4.9786686785484031e+00
39 2.1766057858052013e+00 2.6649010324053237e+00 -5.5648981412947700e+00
40 -8.6278470191715790e+00 -6.1938287573923239e+00 -5.3007116622920181e-01
41 -2.1876614204878253e+00 -4.9711482834141876e+00 2.1023653625053651e+00
42 6.7356638591944380e+00 2.3811284335916083e+00 -2.6936658473425106e+00
43 -4.8166095077507318e+00 9.0096363628437235e+00 -3.5021097239659325e+00
44 4.5856314854049890e-01 -6.4502158514658676e+00 -2.2547412879230078e+00
45 -1.4050653828419732e+00 -2.7741840518241623e-01 5.9455907460299153e+00
46 1.7068092414266975e+00 3.0134108012937979e+00 -7.3633262984858643e+00
47 -8.5993680232001868e+00 3.9482405710461826e-02 -7.4020141073315715e-01
48 -1.4797111864775530e+02 9.0978008317055270e+01 8.8422607641864516e+01
49 1.2110838856400735e+01 1.0210083421020201e+01 -1.7064250014241885e+01
50 1.5673894480674800e+01 4.2472425008217618e+00 2.3874134866435037e-01
51 -1.0755978881146193e+03 9.4584826950708475e+02 -1.2161507490133176e+03
52 1.5238648316158450e+01 -1.4713914840395454e+01 1.0700180235730782e+01
53 1.0832818442423763e+03 -9.4512225095469830e+02 1.2182464933546009e+03
54 -1.2715521453452059e+00 1.0003156043833941e+00 -4.7435368103541586e+00
55 -1.2828505786775382e+01 2.4021051786462557e+01 -1.1214232355595612e+01
56 -7.6186543306062049e+01 2.6038364737554620e+01 4.6298155917897482e+01
57 1.0443716577117205e+02 3.6009056789261075e+02 2.1380802968687473e+02
58 -9.1275826241798370e+00 -1.3245000554098924e+01 -8.7148296348715926e-01
59 1.2840492718836813e+02 -1.0559491971264805e+02 -6.6237876488606588e+01
60 -2.1677737255323635e+00 -1.1892793887650410e+00 -8.1033360541735249e+00
61 -1.2148578617098147e+02 -3.5657463872647236e+02 -2.0401151842776562e+02
62 -2.7445863145191973e+00 1.2181869352734287e+01 6.0139221094528885e+00
63 -7.8700771372090834e+00 -3.1538926482273300e+00 1.0493868436023858e+01
64 2.5378792969187284e+01 3.2092052615283109e+01 -1.6443232536994536e+01
run_vdwl: -1745.867185907274
run_coul: -134.07506093090166
1 4.7247694465962518e+00 3.8636589800094097e+00 1.1621871799061909e+01
2 9.1252156972042364e-02 4.2008105070379278e+00 -1.1980129104691510e+00
3 -1.2601754020721603e+00 9.2524682920861423e+00 7.3929409765368757e+00
4 -2.0726977587259077e+00 2.9077845151833932e+00 -2.7757587746193373e+00
5 -1.0459077110812462e+01 -1.1826792529252666e+01 -1.1083450570197344e+01
6 1.2136543334776942e+01 4.5849120194674704e+00 1.2079301503775022e+00
7 -2.2741379605274665e+00 1.8070366300153946e+00 4.4062289068411875e+00
8 -4.7711645442680819e+00 8.7524859750982991e+00 -4.8962583041964738e+00
9 4.0328267692134583e+00 3.4040816708138846e+00 -1.8995726137368742e+00
10 7.3342970772780651e+00 -1.1625708495684536e+01 1.6641556656164258e+00
11 2.3806275027928057e+00 2.0941569224028513e+00 -1.6519953937630087e+00
12 -6.2260236357885844e+02 8.2842855555945846e+02 -4.9419366949715675e+02
13 3.9701791827302380e+01 2.5107306829177205e+01 -2.6693078280374319e+01
14 -4.9665862218475114e+00 6.7166971146470966e-02 1.0835738495219992e+01
15 -6.5972542674545345e+00 3.2687337906397440e+00 -2.0185347166446563e+00
16 -2.9095101910253285e-01 -7.0769446878294087e+00 -2.6818145980741166e+00
17 -4.6662839821444395e+00 1.8992979157452603e+00 1.3892844154040138e+00
18 3.6232824669878743e+00 -4.3996978419370389e+00 -7.9083611962231020e+00
19 -6.8890131136555359e+00 -1.0830898273912712e+00 1.2072661383321417e+00
20 9.8171454663225788e+00 -9.1035000219673612e-01 2.8436105157483560e+00
21 -6.4267201894937989e+00 -8.2156139232861101e+00 6.3067572266832048e+00
22 -8.5030200202157591e+00 1.4205269387466176e+01 -7.8844799176414240e+00
23 6.9432523243590740e+00 3.3737260499363608e+00 2.2143562043946177e+00
24 1.1167733798586100e+01 -1.0739684016148335e+01 -3.7956397167384850e+00
25 -7.3628620125144280e+01 1.0248504583615836e+02 -4.8418860956676525e+01
26 6.6263639105567663e+01 -6.4097840274469718e+01 -6.3595212300135067e+01
27 7.2057038082439657e+00 7.4557411533732756e+00 2.8696162515374115e+00
28 -7.9329300958368369e+00 1.2335261827158241e+00 -1.2407218911970331e+01
29 -7.9031028432291606e+01 6.4797754884349033e+01 7.2935815813928102e+01
30 -4.3242124829998954e+00 1.7245335206633547e+00 2.5245835369028411e+00
31 -3.2556531400191888e+00 -4.0771319166777253e+00 9.4144942380785412e+00
32 5.3493297950927108e+00 1.2562955135214278e+01 -4.4882693411414065e-01
33 -1.2255425910828938e+01 -2.8689448972424191e+01 1.7093467049309599e+01
34 -3.1777406620351565e+00 -8.9186313551255836e+00 1.0463596621056004e+00
35 -3.1196549028141418e+00 -2.5994883508560860e+00 -9.2177540665900048e-01
36 -5.9576900168680282e+01 -7.4451272068728869e+01 6.1235068709614346e+01
37 6.0157330299178689e+01 8.0285495378129127e+01 -5.7938854098325038e+01
38 -9.6815577392778085e-01 5.8373609530147199e+00 -5.4776558838012690e+00
39 6.5036942143592802e-01 -5.5095290902514762e+00 -8.5045800982305975e-01
40 5.1842254872430082e-01 -1.7683063757591517e+00 -5.3317544095977998e-01
41 -1.6916180625756201e+01 -1.6630288797681160e+01 -3.7719493322749247e+00
42 2.5876740680807822e+00 5.3905848734122879e+00 -1.6450920338360145e+01
43 -3.5279855851679631e+00 9.8495179159335908e+00 -3.2714238772017765e+00
44 6.0889002640825662e-01 3.1410890913935359e+00 1.1032945837449786e+00
45 1.1968416936354009e+01 1.1004153957788535e+01 1.2105055934473947e+01
46 5.9540499922261426e+02 -8.5332666872039056e+02 5.2223674210671663e+02
47 9.3085758478118343e+00 3.6902311757607924e+00 -3.0446774945196609e+00
48 7.5901220066242217e+01 1.3757697268206282e+01 -6.5987345219634847e+01
49 -6.8359329024029364e+01 -3.9765670869480545e+01 5.4460314514030138e+01
50 -4.5653564184344297e-01 -4.6827143163841587e+00 2.2100068640440114e+01
51 -1.0710153088032328e+03 9.5012834910106676e+02 -1.2142374746384971e+03
52 1.3972268755441464e+01 -1.5027397170129637e+01 1.1181740632667148e+01
53 1.0887040837377938e+03 -9.4464616532703872e+02 1.2184820624884132e+03
54 -1.4916779565987026e+01 5.8262752191213307e+00 -1.8719794490582007e+01
55 5.1258703585615493e+01 -1.0313579382868322e+02 3.9556815338722231e+01
56 7.2230447113797958e+00 2.2629738499143413e+00 -3.2153882168284924e+00
57 7.9701310882685341e+00 1.9747232318959853e+01 1.0450527723051438e+01
58 -5.5233659987382948e+00 9.0259835310388870e-01 4.6808879589599970e+00
59 -2.1517171549006143e+01 1.6544754727019644e+01 9.6541292269942787e+00
60 1.2315294601931821e+00 -9.8191775257502307e+00 -1.0095770367374730e+01
61 -1.4403329803680169e+01 -2.5271904205234947e+01 -8.4385127431848996e+00
62 -7.2508809251494712e+00 -4.3099363976614820e+00 3.6806953545058976e+00
63 1.7366944215878561e+01 -2.1678065532175923e+00 -3.9215344710226812e+00
64 2.7331835514835586e+01 2.8927730528687206e+01 -1.7474424636631589e+01
run_vdwl: -2140.050580658317
run_coul: -134.89831768121073
run_stress: ! |-
-4.5872997550021264e+03 -3.9988755027842894e+03 -6.1188498400678936e+03 4.1760773642646591e+03 -5.3312812896810128e+03 4.4665905499115461e+03
-5.4348057135758700e+03 -5.3159515873286382e+03 -6.7051641305464191e+03 5.3862756055333693e+03 -5.9963166859345883e+03 6.0214993931030604e+03
run_forces: ! |2
1 7.9285352946497758e+00 -2.1027294435031489e+00 1.1086825816163195e+01
2 -8.5419611035054892e+00 -1.8545741670495524e+00 7.0079173935199286e+00
3 8.4201002043041484e+00 -2.3236585201563706e+01 -3.5619050765611897e+00
4 5.8244272716577141e+00 2.0336700262604413e+01 -3.7991301665341801e+00
5 -1.0597147068623588e+01 -1.2052626644608324e+01 -8.3616620515931093e+00
6 2.3454925059969298e+00 1.9423650844033959e+00 4.1282352821875090e+00
7 1.1498130683730052e+01 2.7037305516932548e+00 -4.3590618591979624e-01
8 -8.7055079862119680e+00 -8.2080543611619063e+00 3.8502610802694401e+00
9 7.8905280043672006e+00 -9.1748249725951716e+00 -4.3595650830599473e+00
10 -1.1214544607539489e+01 -2.0450291190922241e+00 8.0819592940862215e+00
11 3.1110360595812111e+01 -3.0379448305557091e+01 -8.0086245459647699e+01
12 1.8159177163415904e+01 4.8329111811515117e+00 -1.8760362770452407e+01
13 -8.3967367015110703e+00 -1.1826424656028374e+01 1.0490586963619998e+01
14 -2.7403143901237863e+00 1.5146525228727186e+00 -9.5102191320553011e+00
15 4.6084095676492058e+00 -3.2455057110391716e+00 -2.0076394841639451e+00
16 -5.9971213141367166e+00 -2.1400797487747161e+00 3.8635660450105780e-01
17 -9.2777494354866672e+00 1.1929670372837402e+01 -1.5781112010078653e+00
18 -1.2640996234737609e+00 -1.8706861715670353e+00 -5.0907920358459542e+00
19 4.7710267464485367e+00 -4.0988795718884213e-01 1.9972858346581028e+00
20 -2.2378279249288067e+00 1.3107343895216937e+01 -4.2428037955197233e+00
21 2.5289893324335391e+00 -7.6558240340245325e+00 -1.9069419539704020e+00
22 -1.2728255711560161e+01 2.2449803210470662e+01 -4.4375794812338860e+00
23 2.2126988647283259e+00 -1.9107846463909647e-01 3.6978572940177510e+00
24 -4.9407746531577494e+00 2.3427146491434673e+01 5.0673415871481255e+00
25 2.0382651093146954e+01 -1.5777195827352591e+01 5.8844240401473211e+00
26 -1.2179841258863439e+01 8.7610073871936400e+00 1.1530406923973768e+01
27 2.6588676610969855e+00 5.1896378052166741e+00 5.6387483053555902e+00
28 1.5666189299730341e+01 1.6249815489114361e+00 1.1769883127324578e+01
29 1.6154781213323030e+01 -1.3566033536379088e+01 -9.1472138878710076e+00
30 -8.3472695417921603e+00 7.0562388858473168e+00 6.3337664999339882e+00
31 4.5583726078147091e-02 -6.8079173592470026e+00 8.8191267004599450e-01
32 -3.7970235777293169e+01 2.1379327045743118e+01 6.2117682066075531e+01
33 -1.6353741251628691e+01 -2.2709420959741056e+01 1.6304846015441370e+01
34 9.9972886088735162e+00 2.2643391331974536e+00 -3.4550103310401896e+00
35 7.6455556546064358e+01 -3.7809722812169873e+01 -5.8409719550108399e+01
36 -6.9348410375051216e+01 -8.5084762326803883e+01 8.7066878049116482e+01
37 6.7398486015543440e+01 8.6977799607188729e+01 -7.7261561475806531e+01
38 7.3284526937033183e+00 -7.6017536517090827e+00 -4.4692146119773133e+00
39 2.2039881666125036e+00 2.6712481163988575e+00 -5.4998118879989439e+00
40 -8.1920739102716507e+00 -5.6013665505583772e+00 -6.3906324471289055e-01
41 -9.6914424345020511e-01 -3.6895968531383003e+00 1.1100392846951848e+00
42 6.4344379128480211e+00 2.8919317863081786e+00 -2.0533867506787638e+00
43 -4.8969216688076900e+00 8.9173540382237473e+00 -3.4262725785712571e+00
44 3.9812820211199523e-01 -6.4925904498103826e+00 -2.2914632965466297e+00
45 -1.1701379586492111e+00 -6.8534375447769547e-01 5.2141786324779167e+00
46 1.7082488047076678e+00 2.9689192858013040e+00 -7.3986980836820644e+00
47 -8.4162004926885050e+00 -1.1549845996840341e-01 -6.9949420861709011e-01
48 -1.7140102399704551e+02 1.0954749031941631e+02 1.0215589866028962e+02
49 1.2958089986488774e+01 1.0763913519899868e+01 -1.6248816454541608e+01
50 1.4518431826929746e+01 5.0424709064095214e+00 5.3597634284364259e-02
51 -1.4015254817893208e+04 1.2504756770641432e+04 -1.5967183598544239e+04
52 1.5420056645525479e+01 -1.4998517432233239e+01 1.0961481290482833e+01
53 1.4023035080488416e+04 -1.2504964789519170e+04 1.5967622019398565e+04
54 -1.2443556258779074e+00 1.1935781697358863e+00 -4.7574145417182354e+00
55 -1.3066707495598941e+01 2.3943248755492228e+01 -1.1562546476831322e+01
56 -8.3304235816705628e+01 2.6818208086831827e+01 5.1686194842892000e+01
57 1.5120792842525938e+02 5.1705018806557780e+02 3.0849363663112024e+02
58 -9.2125676676377548e+00 -1.1733522325859413e+01 1.8639253510409390e-02
59 1.5201721677305994e+02 -1.2328671257917705e+02 -8.0267873143196880e+01
60 -2.1359354701455682e+00 -1.0122954230610202e+00 -7.9908260633230173e+00
61 -1.6847550389189405e+02 -5.1524925258885060e+02 -2.9923087564684147e+02
62 -2.7864227120115106e+00 1.2259886975143148e+01 5.9379812236948695e+00
63 -7.9097586153187249e+00 -3.2603634511507762e+00 1.0417329010466950e+01
64 2.5990005859488274e+01 3.2517151166598666e+01 -1.6862446054183820e+01
1 4.7122965463705064e+00 3.9483331580787571e+00 1.1728507090799370e+01
2 6.2540769620511583e-04 4.1975594096677806e+00 -1.2655578889200010e+00
3 -1.2620410935788802e+00 9.3431467733916858e+00 7.4655223258749022e+00
4 -2.0887790508687458e+00 2.8362620671083052e+00 -2.6665635606311202e+00
5 -1.0651156040553424e+01 -1.1870381165237049e+01 -1.1326632090513900e+01
6 1.2091930350658666e+01 4.9253117250743088e+00 1.2569218632223669e+00
7 -2.2466680039472138e+00 1.8285162496432987e+00 4.3642157656708376e+00
8 -4.7964507313245841e+00 8.5783752296325684e+00 -4.7402266739984285e+00
9 4.0921933327781339e+00 3.4422544385699028e+00 -1.9832969878570745e+00
10 7.3261797536542481e+00 -1.1778853782806936e+01 1.6672873842416867e+00
11 2.3262362451564815e+00 2.0302894707130710e+00 -1.5989681942671028e+00
12 -1.8528259976245022e+03 2.6733690793551527e+03 -1.6526784648042444e+03
13 2.9930583975905577e+01 2.0302968352922949e+01 -1.9718778029909050e+01
14 -4.9854520942538061e+00 1.1769885859903972e-01 1.0836441479065130e+01
15 -6.6436375960539298e+00 3.3504771031444984e+00 -2.0103943675023666e+00
16 -8.2826662234260773e-01 -7.7053582800848890e+00 -3.1407715947253125e+00
17 -4.6285853687124510e+00 1.8831032308670554e+00 1.3311155953669225e+00
18 3.6581901150678120e+00 -4.4248478511248424e+00 -7.9261897763565461e+00
19 -6.9775721937943915e+00 -1.0738344822024193e+00 1.3399278795331213e+00
20 9.8542409756528802e+00 -8.3583213038045168e-01 2.9729110796187017e+00
21 -6.3858473251137298e+00 -8.2063778055692573e+00 6.2090714363591655e+00
22 -8.5268899256470476e+00 1.4061397787130135e+01 -7.5959221227930600e+00
23 7.1053916341971437e+00 3.3593780094184811e+00 2.1454774319199919e+00
24 1.1287923752723026e+01 -1.0905349485022397e+01 -3.8423124003720108e+00
25 -7.8878144031307869e+01 1.1144037793856330e+02 -5.2345476930838664e+01
26 7.3556055368024303e+01 -7.0063862890920788e+01 -6.9625728550847427e+01
27 7.2649904809058601e+00 7.4435155311914185e+00 2.8761503856744834e+00
28 -7.8199489945593141e+00 1.6386033020381725e+00 -1.2582811155181844e+01
29 -8.6450147297519081e+01 7.0686776700164231e+01 7.9023187918102835e+01
30 -4.3027444462158391e+00 1.7239104829814220e+00 2.5535790049846039e+00
31 -3.3046500113513737e+00 -4.1204441409767592e+00 9.2623889633129082e+00
32 5.1772833147365276e+00 1.2309544714901417e+01 -5.8472477813045920e-01
33 -1.2649794837326658e+01 -2.9458570076205088e+01 1.7550527253526816e+01
34 -3.1661089296085292e+00 -8.7689043248204115e+00 9.7058611853842613e-01
35 -3.0924038824986666e+00 -2.4489419534639865e+00 -9.3762339249707094e-01
36 -7.1208690920666641e+01 -8.8806793017351026e+01 7.4268951629268130e+01
37 7.1829111376814453e+01 9.4565190186422143e+01 -7.0841892313239953e+01
38 -7.9134637265759067e-01 6.1278427310858481e+00 -5.7340689571314289e+00
39 5.8399235535558125e-01 -5.5767896971177642e+00 -8.8077256557329009e-01
40 3.3991818916474581e-01 -1.5995866916432444e+00 -3.7663653521179696e-01
41 -1.7039926503399553e+01 -1.6822911562713731e+01 -3.6075749889978770e+00
42 2.3106705960272320e+00 5.8882256112675435e+00 -1.6069778691774765e+01
43 -3.4306138444009289e+00 9.9265316014995850e+00 -3.3621117748415665e+00
44 6.4642428588167933e-01 3.0509717980229727e+00 1.1215034345971302e+00
45 1.2022816953757598e+01 1.0808099503778980e+01 1.1641671818467149e+01
46 1.8356549823565331e+03 -2.6920947307474571e+03 1.6746382302672109e+03
47 9.9557809841328968e+00 2.7776852967803318e+00 -3.0712599231891318e+00
48 8.2517343899876806e+01 1.7204532623712410e+01 -7.1545905171861946e+01
49 -7.5263235134837274e+01 -4.2852137596696686e+01 5.9345417348920002e+01
50 -1.7147346578728082e+00 -3.9215028587126870e+00 2.2046250885945959e+01
51 -1.4038860116616930e+04 1.2534673144126318e+04 -1.6004375713307671e+04
52 1.4036626077387192e+01 -1.5290044850585737e+01 1.1474515373499523e+01
53 1.4056802256914347e+04 -1.2530813382537252e+04 1.6008356188605358e+04
54 -1.4964928080874120e+01 5.9952178313910824e+00 -1.8703881764056817e+01
55 5.7171923852348037e+01 -1.1179386471206521e+02 4.3505525189942027e+01
56 7.4502714304370654e+00 2.4147743497707221e+00 -3.3080635176139359e+00
57 8.5439425697364193e+00 2.1826866783942567e+01 1.1739018146458523e+01
58 -5.4915924082367980e+00 1.2005891987344837e+00 4.8298527235029205e+00
59 -2.1827228692900594e+01 1.6434319959036547e+01 9.8746730849843019e+00
60 1.2577696663345079e+00 -9.7685855242912236e+00 -1.0067453390698732e+01
61 -1.4849037643280727e+01 -2.7508969230557387e+01 -9.8137997879873637e+00
62 -7.0890794009016256e+00 -4.5729644069748323e+00 4.0583900180348680e+00
63 1.7531488436729465e+01 -2.0193633823631418e+00 -4.2035828162489199e+00
64 2.8002375179649803e+01 2.9392313693873362e+01 -1.7921068696312695e+01
...

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