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

Merge branch 'upstream' into regression-tests

Browse Source
This commit is contained in:
Trung Nguyen
2024-07-20 16:46:49 -05:00
parent afa44e931b 1ed6e153b6
commit 9676b0cb06
222 changed files with 19077 additions and 3557 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
cmake/CMakeLists.txt
View File

@ -204,7 +204,7 @@ option(BUILD_LAMMPS_GUI "Build and install the LAMMPS GUI" OFF)
# Support using clang-tidy for C++ files with selected options
set(ENABLE_CLANG_TIDY OFF CACHE BOOL "Include clang-tidy processing when compiling")
if(ENABLE_CLANG_TIDY)
set(CMAKE_CXX_CLANG_TIDY "clang-tidy;-checks=-*,performance-trivially-destructible,performance-unnecessary-copy-initialization,performance-unnecessary-value-param,readability-redundant-control-flow,readability-redundant-declaration,readability-redundant-function-ptr-dereference,readability-redundant-member-init,readability-redundant-string-cstr,readability-redundant-string-init,readability-simplify-boolean-expr,readability-static-accessed-through-instance,readability-static-definition-in-anonymous-namespace,modernize-use-override,modernize-use-bool-literals,modernize-use-emplace,modernize-return-braced-init-list,modernize-use-equals-default,modernize-use-equals-delete,modernize-replace-random-shuffle,modernize-deprecated-headers,modernize-use-nullptr,modernize-use-noexcept,modernize-redundant-void-arg;-fix;-header-filter=.*,header-filter=library.h,header-filter=fmt/*.h" CACHE STRING "clang-tidy settings")
set(CMAKE_CXX_CLANG_TIDY "clang-tidy;-checks=-*,performance-trivially-destructible,performance-unnecessary-copy-initialization,performance-unnecessary-value-param,readability-redundant-control-flow,readability-redundant-declaration,readability-redundant-function-ptr-dereference,readability-redundant-member-init,readability-redundant-string-cstr,readability-redundant-string-init,readability-simplify-boolean-expr,readability-static-accessed-through-instance,readability-static-definition-in-anonymous-namespace,readability-qualified-auto,misc-unused-parameters,modernize-deprecated-ios-base-aliases,modernize-loop-convert,modernize-shrink-to-fit,modernize-use-auto,modernize-use-using,modernize-use-override,modernize-use-bool-literals,modernize-use-emplace,modernize-return-braced-init-list,modernize-use-equals-default,modernize-use-equals-delete,modernize-replace-random-shuffle,modernize-deprecated-headers,modernize-use-nullptr,modernize-use-noexcept,modernize-redundant-void-arg;-fix;-header-filter=.*,header-filter=library.h,header-filter=fmt/*.h" CACHE STRING "clang-tidy settings")
else()
unset(CMAKE_CXX_CLANG_TIDY CACHE)
endif()
@ -306,6 +306,7 @@ set(STANDARD_PACKAGES
REACTION
REAXFF
REPLICA
RHEO
RIGID
SCAFACOS
SHOCK
@ -410,6 +411,7 @@ pkg_depends(CG-DNA ASPHERE)
pkg_depends(ELECTRODE KSPACE)
pkg_depends(EXTRA-MOLECULE MOLECULE)
pkg_depends(MESONT MOLECULE)
pkg_depends(RHEO BPM)
# detect if we may enable OpenMP support by default
set(BUILD_OMP_DEFAULT OFF)
@ -550,7 +552,7 @@ else()
endif()
foreach(PKG_WITH_INCL KSPACE PYTHON ML-IAP VORONOI COLVARS ML-HDNNP MDI MOLFILE NETCDF
PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM COMPRESS ML-PACE LEPTON)
PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM COMPRESS ML-PACE LEPTON RHEO)
if(PKG_${PKG_WITH_INCL})
include(Packages/${PKG_WITH_INCL})
endif()
@ -940,6 +942,7 @@ message(STATUS "<<< Compilers and Flags: >>>
-- C++ Compiler: ${CMAKE_CXX_COMPILER}
Type: ${CMAKE_CXX_COMPILER_ID}
Version: ${CMAKE_CXX_COMPILER_VERSION}
C++ Standard: ${CMAKE_CXX_STANDARD}
C++ Flags: ${CMAKE_CXX_FLAGS} ${CMAKE_CXX_FLAGS_${BTYPE}}
Defines: ${DEFINES}")
get_target_property(OPTIONS lammps COMPILE_OPTIONS)

2
cmake/Modules/Packages/ML-QUIP.cmake
View File

@ -27,7 +27,7 @@ if(DOWNLOAD_QUIP)
else()
message(FATAL_ERROR "The ${CMAKE_Fortran_COMPILER_ID} Fortran compiler is not (yet) supported for building QUIP")
endif()
set(temp "${temp}CFLAGS += -fPIC \nCPLUSPLUSFLAGS += -fPIC\nAR_ADD=src\n")
set(temp "${temp}CFLAGS += -fPIC -Wno-return-mismatch \nCPLUSPLUSFLAGS += -fPIC -Wno-return-mismatch\nAR_ADD=src\n")
set(temp "${temp}MATH_LINKOPTS=")
foreach(flag ${BLAS_LIBRARIES})
set(temp "${temp} ${flag}")

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

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

4
cmake/Modules/Testing.cmake
View File

@ -102,9 +102,9 @@ endif()
#######################################
# select code sanitizer options
#######################################
set(ENABLE_SANITIZER "none" CACHE STRING "Select a code sanitizer option (none (default), address, leak, thread, undefined)")
set(ENABLE_SANITIZER "none" CACHE STRING "Select a code sanitizer option (none (default), address, hwaddress, leak, thread, undefined)")
mark_as_advanced(ENABLE_SANITIZER)
set(ENABLE_SANITIZER_VALUES none address leak thread undefined)
set(ENABLE_SANITIZER_VALUES none address hwaddress leak thread undefined)
set_property(CACHE ENABLE_SANITIZER PROPERTY STRINGS ${ENABLE_SANITIZER_VALUES})
validate_option(ENABLE_SANITIZER ENABLE_SANITIZER_VALUES)
string(TOLOWER ${ENABLE_SANITIZER} ENABLE_SANITIZER)

1
cmake/presets/all_off.cmake
View File

@ -82,6 +82,7 @@ set(ALL_PACKAGES
REACTION
REAXFF
REPLICA
RHEO
RIGID
SCAFACOS
SHOCK

1
cmake/presets/all_on.cmake
View File

@ -84,6 +84,7 @@ set(ALL_PACKAGES
REACTION
REAXFF
REPLICA
RHEO
RIGID
SCAFACOS
SHOCK

5
doc/src/Build_development.rst
View File

@ -88,8 +88,8 @@ on recording all commands required to do the compilation.
.. _sanitizer:
Address, Undefined Behavior, and Thread Sanitizer Support (CMake only)
----------------------------------------------------------------------
Address, Leak, Undefined Behavior, and Thread Sanitizer Support (CMake only)
----------------------------------------------------------------------------
Compilers such as GCC and Clang support generating instrumented binaries
which use different sanitizer libraries to detect problems in the code
@ -110,6 +110,7 @@ compilation and linking stages. This is done through setting the
-D ENABLE_SANITIZER=none # no sanitizer active (default)
-D ENABLE_SANITIZER=address # enable address sanitizer / memory leak checker
-D ENABLE_SANITIZER=hwaddress # enable hardware assisted address sanitizer / memory leak checker
-D ENABLE_SANITIZER=leak # enable memory leak checker (only)
-D ENABLE_SANITIZER=undefined # enable undefined behavior sanitizer
-D ENABLE_SANITIZER=thread # enable thread sanitizer

40
doc/src/Build_extras.rst
View File

@ -59,6 +59,7 @@ This is the list of packages that may require additional steps.
* :ref:`POEMS <poems>`
* :ref:`PYTHON <python>`
* :ref:`QMMM <qmmm>`
* :ref:`RHEO <rheo>`
* :ref:`SCAFACOS <scafacos>`
* :ref:`VORONOI <voronoi>`
* :ref:`VTK <vtk>`
@ -1566,10 +1567,11 @@ LAMMPS build.
.. tab:: CMake build
When the ``-D PKG_PLUMED=yes`` flag is included in the cmake
command you must ensure that GSL is installed in locations that
are specified in your environment. There are then two additional
variables that control the manner in which PLUMED is obtained and
linked into LAMMPS.
command you must ensure that `the GNU Scientific Library (GSL)
<https://www.gnu.org/software/gsl/>` is installed in locations
that are accessible in your environment. There are then two
additional variables that control the manner in which PLUMED is
obtained and linked into LAMMPS.
.. code-block:: bash
@ -2040,6 +2042,36 @@ verified to work in February 2020 with Quantum Espresso versions 6.3 to
----------
.. _rheo:
RHEO package
------------
To build with this package you must have the `GNU Scientific Library
(GSL) <https://www.gnu.org/software/gsl/>` installed in locations that
are accessible in your environment. The GSL library should be at least
version 2.7.
.. tabs::
.. tab:: CMake build
If CMake cannot find the GSL library or include files, you can set:
.. code-block:: bash
-D GSL_ROOT_DIR=path # path to root of GSL installation
.. tab:: Traditional make
LAMMPS will try to auto-detect the GSL compiler and linker flags
from the corresponding ``pkg-config`` file (``gsl.pc``), otherwise
you can edit the file ``lib/rheo/Makefile.lammps``
to specify the paths and library names where indicated by comments.
This must be done **before** the package is installed.
----------
.. _scafacos:
SCAFACOS package

8
doc/src/Build_link.rst
View File

@ -45,8 +45,8 @@ executable code from the library is copied into the calling executable.
.. code-block:: bash
mpicc -c -O $(pkgconf liblammps --cflags) caller.c
mpicxx -o caller caller.o -$(pkgconf liblammps --libs)
mpicc -c -O $(pkg-config --cflags liblammps) caller.c
mpicxx -o caller caller.o -$(pkg-config --libs liblammps)
.. tab:: Traditional make
@ -155,8 +155,8 @@ POEMS package installed becomes:
.. code-block:: bash
mpicc -c -O $(pkgconf liblammps --cflags) caller.c
mpicxx -o caller caller.o -$(pkgconf --libs)
mpicc -c -O $(pkg-config --cflags liblammps) caller.c
mpicxx -o caller caller.o -$(pkg-config --libs liblammps)
.. tab:: Traditional make

1
doc/src/Build_package.rst
View File

@ -62,6 +62,7 @@ packages:
* :ref:`POEMS <poems>`
* :ref:`PYTHON <python>`
* :ref:`QMMM <qmmm>`
* :ref:`RHEO <rheo>`
* :ref:`SCAFACOS <scafacos>`
* :ref:`VORONOI <voronoi>`
* :ref:`VTK <vtk>`

1
doc/src/Commands_bond.rst
View File

@ -54,6 +54,7 @@ OPT.
* :doc:`oxdna2/fene <bond_oxdna>`
* :doc:`oxrna2/fene <bond_oxdna>`
* :doc:`quartic (o) <bond_quartic>`
* :doc:`rheo/shell <bond_rheo_shell>`
* :doc:`special <bond_special>`
* :doc:`table (o) <bond_table>`

1
doc/src/Commands_compute.rst
View File

@ -126,6 +126,7 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`reduce <compute_reduce>`
* :doc:`reduce/chunk <compute_reduce_chunk>`
* :doc:`reduce/region <compute_reduce>`
* :doc:`rheo/property/atom <compute_rheo_property_atom>`
* :doc:`rigid/local <compute_rigid_local>`
* :doc:`saed <compute_saed>`
* :doc:`slcsa/atom <compute_slcsa_atom>`

6
doc/src/Commands_fix.rst
View File

@ -28,6 +28,7 @@ OPT.
* :doc:`adapt <fix_adapt>`
* :doc:`adapt/fep <fix_adapt_fep>`
* :doc:`addforce <fix_addforce>`
* :doc:`add/heat <fix_add_heat>`
* :doc:`addtorque <fix_addtorque>`
* :doc:`alchemy <fix_alchemy>`
* :doc:`amoeba/bitorsion <fix_amoeba_bitorsion>`
@ -204,6 +205,11 @@ OPT.
* :doc:`reaxff/species (k) <fix_reaxff_species>`
* :doc:`recenter <fix_recenter>`
* :doc:`restrain <fix_restrain>`
* :doc:`rheo <fix_rheo>`
* :doc:`rheo/oxidation <fix_rheo_oxidation>`
* :doc:`rheo/pressure <fix_rheo_pressure>`
* :doc:`rheo/thermal <fix_rheo_thermal>`
* :doc:`rheo/viscosity <fix_rheo_viscosity>`
* :doc:`rhok <fix_rhok>`
* :doc:`rigid (o) <fix_rigid>`
* :doc:`rigid/meso <fix_rigid_meso>`

6
doc/src/Commands_pair.rst
View File

@ -35,6 +35,10 @@ OPT.
*
*
*
*
*
*
*
* :doc:`adp (ko) <pair_adp>`
* :doc:`agni (o) <pair_agni>`
* :doc:`aip/water/2dm (t) <pair_aip_water_2dm>`
@ -260,6 +264,8 @@ OPT.
* :doc:`rebo (io) <pair_airebo>`
* :doc:`rebomos (o) <pair_rebomos>`
* :doc:`resquared (go) <pair_resquared>`
* :doc:`rheo <pair_rheo>`
* :doc:`rheo/solid <pair_rheo_solid>`
* :doc:`saip/metal (t) <pair_saip_metal>`
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>`
* :doc:`smatb <pair_smatb>`

2
doc/src/Examples.rst
View File

@ -134,6 +134,8 @@ Lowercase directories
+-------------+------------------------------------------------------------------+
| rerun | use of rerun and read_dump commands |
+-------------+------------------------------------------------------------------+
| rheo | RHEO simulations of fluid flows and phase transitions |
+-------------+------------------------------------------------------------------+
| rigid | rigid bodies modeled as independent or coupled |
+-------------+------------------------------------------------------------------+
| shear | sideways shear applied to 2d solid, with and without a void |

1
doc/src/Howto.rst
View File

@ -89,6 +89,7 @@ Packages howto
Howto_drude2
Howto_peri
Howto_manifold
Howto_rheo
Howto_spins
Tutorials howto

18
doc/src/Howto_chunk.rst
View File

@ -1,7 +1,7 @@
Use chunks to calculate system properties
=========================================
In LAMMS, "chunks" are collections of atoms, as defined by the
In LAMMPS, "chunks" are collections of atoms, as defined by the
:doc:`compute chunk/atom <compute_chunk_atom>` command, which assigns
each atom to a chunk ID (or to no chunk at all). The number of chunks
and the assignment of chunk IDs to atoms can be static or change over
@ -148,14 +148,14 @@ Example calculations with chunks
Here are examples using chunk commands to calculate various
properties:
(1) Average velocity in each of 1000 2d spatial bins:
1. Average velocity in each of 1000 2d spatial bins:
.. code-block:: LAMMPS
compute cc1 all chunk/atom bin/2d x 0.0 0.1 y lower 0.01 units reduced
fix 1 all ave/chunk 100 10 1000 cc1 vx vy file tmp.out
(2) Temperature in each spatial bin, after subtracting a flow
2. Temperature in each spatial bin, after subtracting a flow
velocity:
.. code-block:: LAMMPS
@ -164,7 +164,7 @@ velocity:
compute vbias all temp/profile 1 0 0 y 10
fix 1 all ave/chunk 100 10 1000 cc1 temp bias vbias file tmp.out
(3) Center of mass of each molecule:
3. Center of mass of each molecule:
.. code-block:: LAMMPS
@ -172,7 +172,7 @@ velocity:
compute myChunk all com/chunk cc1
fix 1 all ave/time 100 1 100 c_myChunk[*] file tmp.out mode vector
(4) Total force on each molecule and ave/max across all molecules:
4. Total force on each molecule and ave/max across all molecules:
.. code-block:: LAMMPS
@ -183,7 +183,7 @@ velocity:
thermo 1000
thermo_style custom step temp v_xave v_xmax
(5) Histogram of cluster sizes:
5. Histogram of cluster sizes:
.. code-block:: LAMMPS
@ -192,16 +192,16 @@ velocity:
compute size all property/chunk cc1 count
fix 1 all ave/histo 100 1 100 0 20 20 c_size mode vector ave running beyond ignore file tmp.histo
(6) An example for using a per-chunk value to apply per-atom forces to
6. An example for using a per-chunk value to apply per-atom forces to
compress individual polymer chains (molecules) in a mixture, is
explained on the :doc:`compute chunk/spread/atom <compute_chunk_spread_atom>` command doc page.
(7) An example for using one set of per-chunk values for molecule
7. An example for using one set of per-chunk values for molecule
chunks, to create a second set of micelle-scale chunks (clustered
molecules, due to hydrophobicity), is explained on the
:doc:`compute reduce/chunk <compute_reduce_chunk>` command doc page.
(8) An example for using one set of per-chunk values (dipole moment
8. An example for using one set of per-chunk values (dipole moment
vectors) for molecule chunks, spreading the values to each atom in
each chunk, then defining a second set of chunks as spatial bins, and
using the :doc:`fix ave/chunk <fix_ave_chunk>` command to calculate an

7
doc/src/Howto_lammps_gui.rst
View File

@ -571,11 +571,12 @@ General Settings:
size for the text editor and log font of the application can be set.
- *GUI update interval:* Allows to set the time interval between GUI
and data updates during a LAMMPS run in milliseconds. The default is
to update the GUI every 100 milliseconds. This is good for most cases.
For LAMMPS runs that run very fast, however, data may be missed and
to update the GUI every 10 milliseconds. This is good for most cases.
For LAMMPS runs that run *very* fast, however, data may be missed and
through lowering this interval, this can be corrected. However, this
will make the GUI use more resources, which may be a problem on some
computers with slower CPUs. The default value is 100 milliseconds.
computers with slower CPUs and a small number of CPU cores. This
setting may be changed to a value between 1 and 1000 milliseconds.
Accelerators:
^^^^^^^^^^^^^

116
doc/src/Howto_rheo.rst Normal file
View File

@ -0,0 +1,116 @@
Reproducing hydrodynamics and elastic objects (RHEO)
====================================================
The RHEO package is a hybrid implementation of smoothed particle
hydrodynamics (SPH) for fluid flow, which can couple to the :doc:`BPM package
<Howto_bpm>` to model solid elements. RHEO combines these methods to enable
mesh-free modeling of multi-phase material systems. Its SPH solver supports
many advanced options including reproducing kernels, particle shifting, free
surface identification, and solid surface reconstruction. To model fluid-solid
systems, the status of particles can dynamically change between a fluid and
solid state, e.g. during melting/solidification, which determines how they
interact and their physical behavior. The package is designed with modularity
in mind, so one can easily turn various features on/off, adjust physical
details of the system, or develop new capabilities. For instance, the numerics
associated with calculating gradients, reproducing kernels, etc. are separated
into distinctclasses to simplify the development of new integration schemes
which can call these calculations. Additional numerical details can be found in
:ref:`(Palermo) <howto_rheo_palermo>` and
:ref:`(Clemmer) <howto_rheo_clemmer>`.
Note, if you simply want to run a traditional SPH simulation, the :ref:`SPH package
<PKG-SPH>` package is likely better suited for your application. It has fewer advanced
features and therefore benefits from improved performance. The :ref:`MACHDYN
<PKG-MACHDYN>` package for solids may also be relevant for fluid-solid problems.
----------
At the core of the package is :doc:`fix rheo <fix_rheo>` which integrates
particle trajectories and controls many optional features (e.g. the use
of reproducing kernels). In conjunction to fix rheo, one must specify an
instance of :doc:`fix rheo/pressure <fix_rheo_pressure>` and
:doc:`fix rheo/viscosity <fix_rheo_viscosity>` to define a pressure equation
of state and viscosity model, respectively. Optionally, one can model
a heat equation with :doc:`fix rheo/thermal <fix_rheo_thermal>`, which also
allows the user to specify equations for a particle's thermal conductivity,
specific heat, latent heat, and melting temperature. The ordering of these
fixes in an an input script matters. Fix rheo must be defined prior to all
other RHEO fixes.
Typically, RHEO requires atom style rheo. In addition to typical atom
properties like positions and forces, particles store a local density,
viscosity, pressure, and status. If thermal evolution is modeled, one must
use atom style rheo/thermal which also includes a local energy, temperature, and
conductivity. Note that the temperature is always derived from the energy.
This implies the *temperature* attribute of :doc:`the set command <set>` does not
affect particles. Instead, one should use the *sph/e* attribute.
The status variable uses bit-masking to track various properties of a particle
such as its current state of matter (fluid or solid) and its location relative
to a surface. Some of these properties (and others) can be accessed using
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`. The *status*
attribute in :doc:`the set command <set>` only allows control over the first bit
which sets the state of matter, 0 is fluid and 1 is solid.
Fluid interactions, including pressure forces, viscous forces, and heat exchange,
are calculated using :doc:`pair rheo <pair_rheo>`. Unlike typical pair styles,
pair rheo ignores the :doc:`special bond <special_bonds>` settings. Instead,
it determines whether to calculate forces based on the status of particles: e.g.,
hydrodynamic forces are only calculated if a fluid particle is involved.
----------
To model elastic objects, there are currently two mechanisms in RHEO, one designed
for bulk solid bodies and the other for thin shells. Both mechanisms rely on
introducing bonded forces between particles and therefore require a hybrid of atom
style bond and rheo (or rheo/thermal).
To create an elastic solid body, one has to (a) change the status of constituent
particles to solid (e.g. with the :doc:`set <set>` command), (b) create bpm
bonds between the particles (see the :doc:`bpm howto <Howto_bpm>` page for
more details), and (c) use :doc:`pair rheo/solid <pair_rheo_solid>` to
apply repulsive contact forces between distinct solid bodies. Akin to pair rheo,
pair rheo/solid considers a particles fluid/solid phase to determine whether to
apply forces. However, unlike pair rheo, pair rheo/solid does obey special bond
settings such that contact forces do not have to be calculated between two bonded
solid particles in the same elastic body.
In systems with thermal evolution, fix rheo/thermal can optionally set a
melting/solidification temperature allowing particles to dynamically swap their
state between fluid and solid when the temperature exceeds or drops below the
critical temperature, respectively. Using the *react* option, one can specify a maximum
bond length and a bond type. Then, when solidifying, particles will search their
local neighbors and automatically create bonds with any neighboring solid particles
in range. For BPM bond styles, bonds will then use the immediate position of the two
particles to calculate a reference state. When melting, particles will delete any
bonds of the specified type when reverting to a fluid state. Special bonds are updated
as bonds are created/broken.
The other option for elastic objects is an elastic shell that is nominally much
thinner than a particle diameter, e.g. a oxide skin which gradually forms over time
on the surface of a fluid. Currently, this is implemented using
:doc:`fix rheo/oxidation <fix_rheo_oxidation>` and bond style
:doc:`rheo/shell <bond_rheo_shell>`. Essentially, fix rheo/oxidation creates candidate
bonds of a specified type between surface fluid particles within a specified distance.
a newly created rheo/shell bond will then start a timer. While the timer is counting
down, the bond will delete itself if particles move too far apart or move away from the
surface. However, if the timer reaches a user-defined threshold, then the bond will
activate and apply additional forces to the fluid particles. Bond style rheo/shell
then operates very similarly to a BPM bond style, storing a reference length and
breaking if stretched too far. Unlike the above method, this option does not remove
the underlying fluid interactions (although particle shifting is turned off) and does
not modify special bond settings of particles.
While these two options are not expected to be appropriate for every system,
either framework can be modified to create more suitable models (e.g. by changing the
criteria for creating/deleting a bond or altering force calculations).
----------
.. _howto_rheo_palermo:
**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, in preparation.
.. _howto_rheo_clemmer:
**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

71
doc/src/Packages_details.rst
View File

@ -8,12 +8,12 @@ info on how to download or build any extra library it requires. It also
gives links to documentation, example scripts, and pictures/movies (if
available) that illustrate use of the package.
The majority of packages can be included in a LAMMPS build with a
single setting (``-D PKG_<NAME>=on`` for CMake) or command
(``make yes-<name>`` for make). See the :doc:`Build package <Build_package>`
page for more info. A few packages may require additional steps;
this is indicated in the descriptions below. The :doc:`Build extras <Build_extras>`
page gives those details.
The majority of packages can be included in a LAMMPS build with a single
setting (``-D PKG_<NAME>=on`` for CMake) or command (``make yes-<name>``
for make). See the :doc:`Build package <Build_package>` page for more
info. A few packages may require additional steps; this is indicated in
the descriptions below. The :doc:`Build extras <Build_extras>` page
gives those details.
.. note::
@ -103,6 +103,7 @@ page gives those details.
* :ref:`QEQ <PKG-QEQ>`
* :ref:`QMMM <PKG-QMMM>`
* :ref:`QTB <PKG-QTB>`
* :ref:`RHEO <PKG-RHEO>`
* :ref:`REACTION <PKG-REACTION>`
* :ref:`REAXFF <PKG-REAXFF>`
* :ref:`REPLICA <PKG-REPLICA>`
@ -1323,18 +1324,19 @@ KSPACE package
**Contents:**
A variety of long-range Coulombic solvers, as well as pair styles
which compute the corresponding short-range pairwise Coulombic
interactions. These include Ewald, particle-particle particle-mesh
(PPPM), and multilevel summation method (MSM) solvers.
A variety of long-range Coulombic solvers, as well as pair styles which
compute the corresponding short-range pairwise Coulombic interactions.
These include Ewald, particle-particle particle-mesh (PPPM), and
multilevel summation method (MSM) solvers.
**Install:**
Building with this package requires a 1d FFT library be present on
your system for use by the PPPM solvers. This can be the KISS FFT
library provided with LAMMPS, third party libraries like FFTW, or a
vendor-supplied FFT library. See the :doc:`Build settings <Build_settings>` page for details on how to select
different FFT options for your LAMPMS build.
Building with this package requires a 1d FFT library be present on your
system for use by the PPPM solvers. This can be the KISS FFT library
provided with LAMMPS, third party libraries like FFTW, or a
vendor-supplied FFT library. See the :doc:`Build settings
<Build_settings>` page for details on how to select different FFT
options for your LAMMPS build.
**Supporting info:**
@ -2621,6 +2623,45 @@ another set.
----------
.. _PKG-RHEO:
RHEO package
------------
**Contents:**
Pair styles, bond styles, fixes, and computes for reproducing
hydrodynamics and elastic objects. See the :doc:`Howto rheo
<Howto_rheo>` page for an overview.
**Install:**
This package has :ref:`specific installation instructions <rheo>` on the :doc:`Build extras <Build_extras>` page.
**Authors:** Joel T. Clemmer (Sandia National Labs),
Thomas C. O'Connor (Carnegie Mellon University)
.. versionadded:: TBD
**Supporting info:**
* src/RHEO filenames -> commands
* :doc:`Howto_rheo <Howto_rheo>`
* :doc:`atom_style rheo <atom_style>`
* :doc:`atom_style rheo/thermal <atom_style>`
* :doc:`bond_style rheo/shell <bond_rheo_shell>`
* :doc:`compute rheo/property/atom <compute_rheo_property_atom>`
* :doc:`fix rheo <fix_rheo>`
* :doc:`fix rheo/oxidation <fix_rheo_oxidation>`
* :doc:`fix rheo/pressure <fix_rheo_pressure>`
* :doc:`fix rheo/thermal <fix_rheo_thermal>`
* :doc:`fix rheo/viscosity <fix_rheo_viscosity>`
* :doc:`pair_style rheo <pair_rheo>`
* :doc:`pair_style rheo/solid <pair_rheo_solid>`
* examples/rheo
----------
.. _PKG-RIGID:
RIGID package

5
doc/src/Packages_list.rst
View File

@ -413,6 +413,11 @@ whether an extra library is needed to build and use the package:
- :doc:`fix qtb <fix_qtb>` :doc:`fix qbmsst <fix_qbmsst>`
- qtb
- no
* - :ref:`RHEO <PKG-RHEO>`
- reproducing hydrodynamics and elastic objects
- :doc:`Howto rheo <Howto_rheo>`
- rheo
- no
* - :ref:`REACTION <PKG-REACTION>`
- chemical reactions in classical MD
- :doc:`fix bond/react <fix_bond_react>`

6
doc/src/Tools.rst
View File

@ -1329,7 +1329,7 @@ for Tcl with:
.. code-block:: bash
swig -tcl -module tcllammps lammps.i
gcc -fPIC -shared $(pkgconf --cflags tcl) -o tcllammps.so \
gcc -fPIC -shared $(pkg-config tcl --cflags) -o tcllammps.so \
lammps_wrap.c -L ../src/ -llammps
tclsh
@ -1340,8 +1340,8 @@ functions included with:
swig -tcl -module tcllmps lammps_shell.i
gcc -o tcllmpsh lammps_wrap.c -Xlinker -export-dynamic \
-DHAVE_CONFIG_H $(pkgconf --cflags tcl) \
$(pkgconf --libs tcl) -L ../src -llammps
-DHAVE_CONFIG_H $(pkg-config tcl --cflags) \
$(pkg-config tcl --libs) -L ../src -llammps
In both cases it is assumed that the LAMMPS library was compiled
as a shared library in the ``src`` folder. Otherwise the last

2
doc/src/atom_modify.rst
View File

@ -190,7 +190,7 @@ Default
By default, *id* is yes. By default, atomic systems (no bond topology
info) do not use a map. For molecular systems (with bond topology
info), the default is to use a map of either *array* or *hash* style
depending on the size of the sustem, as explained above for the *map
depending on the size of the system, as explained above for the *map
yes* keyword/value option. By default, a *first* group is not
defined. By default, sorting is enabled with a frequency of 1000 and
a binsize of 0.0, which means the neighbor cutoff will be used to set

8
doc/src/atom_style.rst
View File

@ -189,6 +189,14 @@ the Additional Information section below.
- *atomic* + molecule, radius, rmass + "smd data"
- :ref:`MACHDYN <PKG-MACHDYN>`
- Smooth Mach Dynamics models
* - *rheo*
- *atomic* + rho, status
- :ref:`RHEO <PKG-RHEO>`
- solid and fluid RHEO particles
* - *rheo/thermal*
- *atomic* + rho, status, energy, temperature
- :ref:`RHEO <PKG-RHEO>`
- RHEO particles with temperature
* - *sph*
- *atomic* + "sph data"
- :ref:`SPH <PKG-SPH>`

188
doc/src/bond_rheo_shell.rst Normal file
View File

@ -0,0 +1,188 @@
.. index:: bond_style rheo/shell
bond_style rheo/shell command
=============================
Syntax
""""""
.. code-block:: LAMMPS
bond_style rheo/shell keyword value attribute1 attribute2 ...
* required keyword = *t/form*
* optional keyword = *store/local*
.. parsed-literal::
*t/form* value = formation time for a bond (time units)
*store/local* values = fix_ID N attributes ...
* fix_ID = ID of associated internal fix to store data
* N = prepare data for output every this many timesteps
* attributes = zero or more of the below attributes may be appended
*id1, id2* = IDs of 2 atoms in the bond
*time* = the timestep the bond broke
*x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units)
*x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units)
Examples
""""""""
.. code-block:: LAMMPS
bond_style rheo/shell t/form 10.0
bond_coeff 1 1.0 0.05 0.1
Description
"""""""""""
.. versionadded:: TBD
The *rheo/shell* bond style is designed to work with
:doc:`fix rheo/oxidation <fix_rheo_oxidation>` which creates candidate
bonds between eligible surface or near-surface particles. When a bond
is first created, it computes no forces and starts a timer. Forces are
not computed until the timer reaches the specified bond formation time,
*t/form*, and the bond is enabled and applies forces. If the two particles
move outside of the maximum bond distance or move into the bulk before
the timer reaches *t/form*, the bond automatically deletes itself. This
deletion is not recorded as a broken bond in the optional *store/local* fix.
Before bonds are enabled, they are still treated as regular bonds by
all other parts of LAMMPS. This means they are written to data files
and counted in computes such as :doc:`nbond/atom <compute_nbond_atom>`.
To only count enabled bonds, use the *nbond/shell* attribute in
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`.
When enabled, the bond then computes forces based on deviations from
the initial reference state of the two atoms much like a BPM style
bond (as further discussed in the :doc:`BPM howto page <Howto_bpm>`).
The reference state is stored by each bond when it is first enabled.
Data is then preserved across run commands and is written to
:doc:`binary restart files <restart>` such that restarting the system
will not reset the reference state of a bond or the timer.
This bond style is based on a model described in
:ref:`(Clemmer) <rheo_clemmer>`. The force has a magnitude of
.. math::
F = 2 k (r - r_0) + \frac{2 k}{r_0^2 \epsilon_c^2} (r - r_0)^3
where :math:`k` is a stiffness, :math:`r` is the current distance
and :math:`r_0` is the initial distance between the two particles, and
:math:`\epsilon_c` is maximum strain beyond which a bond breaks. This
is done by setting the bond type to 0 such that forces are no longer
computed.
A damping force proportional to the difference in the normal velocity
of particles is also applied to bonded particles:
.. math::
F_D = - \gamma w (\hat{r} \bullet \vec{v})
where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the
displacement normal vector, and :math:`\vec{v}` is the velocity difference
between the two particles.
The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data
<read_data>` or :doc:`read_restart <read_restart>` commands:
* :math:`k` (force/distance units)
* :math:`\epsilon_c` (unit less)
* :math:`\gamma` (force/velocity units)
Unlike other BPM-style bonds, this bond style does not update special
bond settings when bonds are created or deleted. This bond style also
does not enforce specific :doc:`special_bonds <special_bonds>` settings.
This behavior is purposeful such :doc:`RHEO pair <pair_rheo>` forces
and heat flows are still calculated.
If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed
and transferred to an internal fix labeled *fix_ID*. This allows the
local data to be accessed by other LAMMPS commands. Following this optional
keyword, a list of one or more attributes is specified. These include the
IDs of the two atoms in the bond. The other attributes for the two atoms
include the timestep during which the bond broke and the current/initial
center of mass position of the two atoms.
Data is continuously accumulated over intervals of *N*
timesteps. At the end of each interval, all of the saved accumulated
data is deleted to make room for new data. Individual datum may
therefore persist anywhere between *1* to *N* timesteps depending on
when they are saved. This data can be accessed using the *fix_ID* and a
:doc:`dump local <dump>` command. To ensure all data is output,
the dump frequency should correspond to the same interval of *N*
timesteps. A dump frequency of an integer multiple of *N* can be used
to regularly output a sample of the accumulated data.
Note that when unbroken bonds are dumped to a file via the
:doc:`dump local <dump>` command, bonds with type 0 (broken bonds)
are not included.
The :doc:`delete_bonds <delete_bonds>` command can also be used to
query the status of broken bonds or permanently delete them, e.g.:
.. code-block:: LAMMPS
delete_bonds all stats
delete_bonds all bond 0 remove
----------
Restart and other info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This bond style writes the reference state of each bond to
:doc:`binary restart files <restart>`. Loading a restart
file will properly restore bonds. However, the reference state is NOT
written to data files. Therefore reading a data file will not
restore bonds and will cause their reference states to be redefined.
If the *store/local* option is used, an internal fix will calculate
a local vector or local array depending on the number of input values.
The length of the vector or number of rows in the array is the number
of recorded, broken bonds. If a single input is specified, a local
vector is produced. If two or more inputs are specified, a local array
is produced where the number of columns = the number of inputs. The
vector or array can be accessed by any command that uses local values
from a compute as input. See the :doc:`Howto output <Howto_output>` page
for an overview of LAMMPS output options.
The vector or array will be floating point values that correspond to
the specified attribute.
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. The single() function also calculates two
extra bond quantities, the initial distance :math:`r_0` and a time.
These extra quantities can be accessed by the
:doc:`compute bond/local <compute_bond_local>` command as *b1* and *b2*\ .
Restrictions
""""""""""""
This bond style is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`bond_coeff <bond_coeff>`, :doc:`fix rheo/oxidation <fix_rheo_oxidation>`
Default
"""""""
NA
----------
.. _rheo_clemmer:
**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

1
doc/src/bond_style.rst
View File

@ -105,6 +105,7 @@ accelerated styles exist.
* :doc:`oxdna2/fene <bond_oxdna>` - same as oxdna but used with different pair styles
* :doc:`oxrna2/fene <bond_oxdna>` - modified FENE bond suitable for RNA modeling
* :doc:`quartic <bond_quartic>` - breakable quartic bond
* :doc:`rheo/shell <bond_rheo_shell>` - shell bond for oxidation modeling in RHEO
* :doc:`special <bond_special>` - enable special bond exclusions for 1-5 pairs and beyond
* :doc:`table <bond_table>` - tabulated by bond length

1
doc/src/compute.rst
View File

@ -290,6 +290,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
* :doc:`reduce <compute_reduce>` - combine per-atom quantities into a single global value
* :doc:`reduce/chunk <compute_reduce_chunk>` - reduce per-atom quantities within each chunk
* :doc:`reduce/region <compute_reduce>` - same as compute reduce, within a region
* :doc:`rheo/property/atom <compute_rheo_property_atom>` - convert atom attributes in RHEO package to per-atom vectors/arrays
* :doc:`rigid/local <compute_rigid_local>` - extract rigid body attributes
* :doc:`saed <compute_saed>` - electron diffraction intensity on a mesh of reciprocal lattice nodes
* :doc:`slcsa/atom <compute_slcsa_atom>` - perform Supervised Learning Crystal Structure Analysis (SL-CSA)

13
doc/src/compute_nbond_atom.rst
View File

@ -8,10 +8,17 @@ Syntax
.. code-block:: LAMMPS
compute ID group-ID nbond/atom
compute ID group-ID nbond/atom keyword value
* ID, group-ID are documented in :doc:`compute <compute>` command
* nbond/atom = style name of this compute command
* zero or more keyword/value pairs may be appended
* keyword = *bond/type*
.. parsed-literal::
*bond/type* value = *btype*
*btype* = bond type included in count
Examples
""""""""
@ -19,6 +26,7 @@ Examples
.. code-block:: LAMMPS
compute 1 all nbond/atom
compute 1 all nbond/atom bond/type 2
Description
"""""""""""
@ -31,6 +39,9 @@ the :doc:`Howto broken bonds <Howto_bpm>` page for more information.
The number of bonds will be zero for atoms not in the specified
compute group. This compute does not depend on Newton bond settings.
If the keyword *bond/type* is specified, only bonds of *btype* are
counted.
Output info
"""""""""""

143
doc/src/compute_rheo_property_atom.rst Normal file
View File

@ -0,0 +1,143 @@
.. index:: compute rheo/property/atom
compute rheo/property/atom command
==================================
Syntax
""""""
.. code-block:: LAMMPS
compute ID group-ID rheo/property/atom input1 input2 ...
* ID, group-ID are documented in :doc:`compute <compute>` command
* rheo/property/atom = style name of this compute command
* input = one or more atom attributes
.. parsed-literal::
possible attributes = phase, surface, surface/r,
surface/divr, surface/n/a, coordination,
shift/v/a, energy, temperature, heatflow,
conductivity, cv, viscosity, pressure, rho,
grad/v/ab, stress/v/ab, stress/t/ab, nbond/shell
.. parsed-literal::
*phase* = atom phase state
*surface* = atom surface status
*surface/r* = atom distance from the surface
*surface/divr* = divergence of position at atom position
*surface/n/a* = a-component of surface normal vector
*coordination* = coordination number
*shift/v/a* = a-component of atom shifting velocity
*energy* = atom energy
*temperature* = atom temperature
*heatflow* = atom heat flow
*conductivity* = atom conductivity
*cv* = atom specific heat
*viscosity* = atom viscosity
*pressure* = atom pressure
*rho* = atom density
*grad/v/ab* = ab-component of atom velocity gradient tensor
*stress/v/ab* = ab-component of atom viscous stress tensor
*stress/t/ab* = ab-component of atom total stress tensor (pressure and viscous)
*nbond/shell* = number of oxide bonds
Examples
""""""""
.. code-block:: LAMMPS
compute 1 all rheo/property/atom phase surface/r surface/n/* pressure
compute 2 all rheo/property/atom shift/v/x grad/v/xx stress/v/*
Description
"""""""""""
.. versionadded:: TBD
Define a computation that stores atom attributes specific to the RHEO
package for each atom in the group. This is useful so that the values
can be used by other :doc:`output commands <Howto_output>` that take
computes as inputs. See for example, the
:doc:`compute reduce <compute_reduce>`,
:doc:`fix ave/atom <fix_ave_atom>`,
:doc:`fix ave/histo <fix_ave_histo>`,
:doc:`fix ave/chunk <fix_ave_chunk>`, and
:doc:`atom-style variable <variable>` commands.
For vector attributes, e.g. *shift/v/*:math:`\alpha`, one must specify
:math:`\alpha` as the *x*, *y*, or *z* component, e.g. *shift/v/x*.
Alternatively, a wild card \* will include all components, *x* and *y* in
2D or *x*, *y*, and *z* in 3D.
For tensor attributes, e.g. *grad/v/*:math:`\alpha \beta`, one must specify
both :math:`\alpha` and :math:`\beta` as *x*, *y*, or *z*, e.g. *grad/v/xy*.
Alternatively, a wild card \* will include all components. In 2D, this
includes *xx*, *xy*, *yx*, and *yy*. In 3D, this includes *xx*, *xy*, *xz*,
*yx*, *yy*, *yz*, *zx*, *zy*, and *zz*.
Many properties require their respective fixes, listed below in related
commands, be defined. For instance, the *viscosity* attribute is the
viscosity of a particle calculated by
:doc:`fix rheo/viscous <fix_rheo_viscosity>`. The meaning of less obvious
properties is described below.
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
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
have a value of 0. The *surface/n/*:math:`\alpha` properties are the
components of the surface normal vector.
The *shift/v/*:math:`\alpha` properties are the components of the shifting
velocity produced by the *shift* option of :doc:`fix rheo <fix_rheo>`.
The *nbond/shell* property is the number of shell bonds that have been
activated from :doc:`bond style rheo/shell <bond_rheo_shell>`.
The values are stored in a per-atom vector or array as discussed
below. Zeroes are stored for atoms not in the specified group or for
quantities that are not defined for a particular particle in the group
Output info
"""""""""""
This compute calculates a per-atom vector or per-atom array depending
on the number of input values. Generally, if a single input is specified,
a per-atom vector is produced. If two or more inputs are specified, a
per-atom array is produced where the number of columns = the number of
inputs. However, if a wild card \* is used for a vector or tensor, then
the number of inputs is considered to be incremented by the dimension or
the dimension squared, respectively. The vector or array can be accessed
by any command that uses per-atom values from a compute as input. See the
:doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
options.
The vector or array values will be in whatever :doc:`units <units>` the
corresponding attribute is in (e.g., density units for *rho*).
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`dump custom <dump>`, :doc:`compute reduce <compute_reduce>`,
:doc:`fix ave/atom <fix_ave_atom>`, :doc:`fix ave/chunk <fix_ave_chunk>`,
:doc:`fix rheo/viscosity <fix_rheo_viscosity>`,
:doc:`fix rheo/pressure <fix_rheo_pressure>`,
:doc:`fix rheo/thermal <fix_rheo_thermal>`,
:doc:`fix rheo/oxdiation <fix_rheo_oxidation>`,
:doc:`fix rheo <fix_rheo>`
Default
"""""""
none

6
doc/src/fix.rst
View File

@ -193,6 +193,7 @@ accelerated styles exist.
* :doc:`adapt <fix_adapt>` - change a simulation parameter over time
* :doc:`adapt/fep <fix_adapt_fep>` - enhanced version of fix adapt
* :doc:`addforce <fix_addforce>` - add a force to each atom
* :doc:`add/heat <fix_add_heat>` - add a heat flux to each atom
* :doc:`addtorque <fix_addtorque>` - add a torque to a group of atoms
* :doc:`alchemy <fix_alchemy>` - perform an "alchemical transformation" between two partitions
* :doc:`amoeba/bitorsion <fix_amoeba_bitorsion>` - torsion/torsion terms in AMOEBA force field
@ -369,6 +370,11 @@ accelerated styles exist.
* :doc:`reaxff/species <fix_reaxff_species>` - write out ReaxFF molecule information
* :doc:`recenter <fix_recenter>` - constrain the center-of-mass position of a group of atoms
* :doc:`restrain <fix_restrain>` - constrain a bond, angle, dihedral
* :doc:`rheo <fix_rheo>` - integrator for the RHEO package
* :doc:`rheo/thermal <fix_rheo_thermal>` - thermal integrator for the RHEO package
* :doc:`rheo/oxidation <fix_rheo_oxidation>` - create oxidation bonds for the RHEO package
* :doc:`rheo/pressure <fix_rheo_pressure>` - pressure calculation for the RHEO package
* :doc:`rheo/viscosity <fix_rheo_pressure>` - viscosity calculation for the RHEO package
* :doc:`rhok <fix_rhok>` - add bias potential for long-range ordered systems
* :doc:`rigid <fix_rigid>` - constrain one or more clusters of atoms to move as a rigid body with NVE integration
* :doc:`rigid/meso <fix_rigid_meso>` - constrain clusters of mesoscopic SPH/SDPD particles to move as a rigid body

111
doc/src/fix_add_heat.rst Normal file
View File

@ -0,0 +1,111 @@
.. index:: fix add/heat
fix add/heat command
====================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID add/heat style args keyword values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* add/heat = style name of this fix command
* style = *constant* or *linear* or *quartic*
.. parsed-literal::
*constant* args = *rate*
*rate* = rate of heat flow (energy/time units)
*linear* args = :math:`T_{target}` *k*
:math:`T_{target}` = target temperature (temperature units)
*k* = prefactor (energy/(time*temperature) units)
*quartic* args = :math:`T_{target}` *k*
:math:`T_{target}` = target temperature (temperature units)
*k* = prefactor (energy/(time*temperature^4) units)
* zero or more keyword/value pairs may be appended to args
* keyword = *overwrite*
.. parsed-literal::
*overwrite* value = *yes* or *no*
*yes* = sets current heat flow of particle
*no* = adds to current heat flow of particle
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all add/heat constant v_heat
fix 1 all add/heat linear 10.0 1.0 overwrite yes
Description
"""""""""""
This fix adds heat to particles with the temperature attribute every timestep.
Note that this is an internal temperature of a particle intended for use with
non-atomistic models like the discrete element method.
For the *constant* style, heat is added at the specified rate. For the *linear* style,
heat is added at a rate of :math:`k (T_{target} - T)` where :math:`k` is the
specified prefactor, :math:`T_{target}` is the specified target temperature, and
:math:`T` is the temperature of the atom. This may be more representative of a
conductive process. For the *quartic* style, heat is added at a rate of
:math:`k (T_{target}^4 - T^4)`, akin to radiative heat transfer.
The rate or temperature can be can be specified as an equal-style or atom-style
:doc:`variable <variable>`. If the value is a variable, it should be
specified as v_name, where name is the variable name. In this case, the
variable will be evaluated each time step, and its value will be used to
determine the rate of heat added.
Equal-style variables can specify formulas with various mathematical
functions and include :doc:`thermo_style <thermo_style>` command
keywords for the simulation box parameters, time step, and elapsed time
to specify time-dependent heating.
Atom-style variables can specify the same formulas as equal-style
variables but can also include per-atom values, such as atom
coordinates to specify spatially-dependent heating.
If the *overwrite* keyword is set to *yes*, this fix will set the total
heat flow on a particle every timestep, overwriting contributions from pair
styles or other fixes. If *overwrite* is *no*, this fix will add heat on
top of other contributions.
----------
Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
No global or per-atom quantities are stored by this fix for access by various
:doc:`output commands <Howto_output>`. No parameter of this fix can be used
with the *start/stop* keywords of the :doc:`run <run>` command. This fix is
not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This pair style is part of the GRANULAR package. It is
only enabled if LAMMPS was built with that package.
See the :doc:`Build package <Build_package>` page for more info.
This fix requires that atoms store temperature and heat flow
as defined by the :doc:`fix property/atom <fix_property_atom>` command.
Related commands
""""""""""""""""
:doc:`fix heat/flow <fix_heat_flow>`,
:doc:`fix property/atom <fix_property_atom>`,
:doc:`fix rheo/thermal <fix_rheo_thermal>`
Default
"""""""
The default for the *overwrite* keyword is *no*

66
doc/src/fix_deform.rst
View File

@ -4,9 +4,6 @@
fix deform command
==================
:doc:`fix deform/pressure <fix_deform_pressure>` command
========================================================
Accelerator Variants: *deform/kk*
Syntax
@ -14,12 +11,11 @@ Syntax
.. code-block:: LAMMPS
fix ID group-ID fix_style N parameter style args ... keyword value ...
fix ID group-ID deform N parameter style args ... keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* fix_style = *deform* or *deform/pressure*
* N = perform box deformation every this many timesteps
* one or more parameter/style/args sequences of arguments may be appended
* one or more parameter/args sequences may be appended
.. parsed-literal::
@ -46,12 +42,6 @@ Syntax
*variable* values = v_name1 v_name2
v_name1 = variable with name1 for box length change as function of time
v_name2 = variable with name2 for change rate as function of time
*pressure* values = target gain (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
*pressure/mean* values = target gain (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
*xy*, *xz*, *yz* args = style value
style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable*
@ -64,8 +54,6 @@ Syntax
effectively an engineering shear strain rate
*erate* value = R
R = engineering shear strain rate (1/time units)
*erate/rescale* value = R (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
R = engineering shear strain rate (1/time units)
*trate* value = R
R = true shear strain rate (1/time units)
*wiggle* values = A Tp
@ -74,9 +62,6 @@ Syntax
*variable* values = v_name1 v_name2
v_name1 = variable with name1 for tilt change as function of time
v_name2 = variable with name2 for change rate as function of time
*pressure* values = target gain (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
* zero or more keyword/value pairs may be appended
* keyword = *remap* or *flip* or *units* or *couple* or *vol/balance/p* or *max/rate* or *normalize/pressure*
@ -92,15 +77,6 @@ Syntax
*units* value = *lattice* or *box*
lattice = distances are defined in lattice units
box = distances are defined in simulation box units
*couple* value = *none* or *xyz* or *xy* or *yz* or *xz* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
couple pressure values of various dimensions
*vol/balance/p* value = *yes* or *no* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
Modifies the behavior of the *volume* option to try and balance pressures
*max/rate* value = *rate* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
rate = maximum strain rate for pressure control
*normalize/pressure* value = *yes* or *no* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
Modifies pressure controls such that the deviation in pressure is normalized by the target pressure
Examples
""""""""
@ -112,8 +88,6 @@ Examples
fix 1 all deform 1 xy erate 0.001 remap v
fix 1 all deform 10 y delta -0.5 0.5 xz vel 1.0
See examples for :doc:`fix deform/pressure <fix_deform_pressure>` on its doc page
Description
"""""""""""
@ -123,17 +97,13 @@ run. Orthogonal simulation boxes have 3 adjustable parameters
adjustable parameters (x,y,z,xy,xz,yz). Any or all of them can be
adjusted independently and simultaneously.
The fix deform command allows use of all the arguments listed above,
except those flagged as available ONLY for the :doc:`fix
deform/pressure <fix_deform_pressure>` command, which are
pressure-based controls. The fix deform/pressure command allows use
of all the arguments listed above.
The rest of this doc page explains the options common to both
commands. The :doc:`fix deform/pressure <fix_deform_pressure>` doc
page explains the options available ONLY with the fix deform/pressure
command. Note that a simulation can define only a single deformation
command: fix deform or fix deform/pressure.
The :doc:`fix deform/pressure <fix_deform_pressure>` command extends
this command with additional keywords and arguments. The rest of this
page explains the options common to both commands. The :doc:`fix
deform/pressure <fix_deform_pressure>` page explains the options
available ONLY with the fix deform/pressure command. Note that a
simulation can define only a single deformation command: fix deform or
fix deform/pressure.
Both these fixes can be used to perform non-equilibrium MD (NEMD)
simulations of a continuously strained system. See the :doc:`fix
@ -143,6 +113,24 @@ simulation of a continuously extended system (extensional flow) can be
modeled using the :ref:`UEF package <PKG-UEF>` and its :doc:`fix
commands <fix_nh_uef>`.
.. admonition:: Inconsistent trajectories due to image flags
:class: warning
When running long simulations while shearing the box or using a high
shearing rate, it is possible that the image flags used for storing
unwrapped atom positions will "wrap around". When LAMMPS is compiled
with the default settings, case image flags are limited to a range of
:math:`-512 \le i \le 511`, which will overflow when atoms starting
at zero image flag value have passed through a periodic box dimension
more than 512 times.
Changing the :ref:`size of LAMMPS integer types <size>` to the
"bigbig" setting can make this overflow much less likely, since it
increases the image flag value range to :math:`- 1,048,576 \le i \le
1\,048\,575`
----------
For the *x*, *y*, *z* parameters, the associated dimension cannot be
shrink-wrapped. For the *xy*, *yz*, *xz* parameters, the associated
second dimension cannot be shrink-wrapped. Dimensions not varied by

71
doc/src/fix_deform_pressure.rst
View File

@ -13,29 +13,66 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* deform/pressure = style name of this fix command
* N = perform box deformation every this many timesteps
* one or more parameter/arg sequences may be appended
* one or more parameter/args sequences may be appended
.. parsed-literal::
parameter = *x* or *y* or *z* or *xy* or *xz* or *yz* or *box*
*x*, *y*, *z* args = style value(s)
style = *final* or *delta* or *scale* or *vel* or *erate* or *trate* or *volume* or *wiggle* or *variable* or *pressure* or *pressure/mean*
*final* values = lo hi
lo hi = box boundaries at end of run (distance units)
*delta* values = dlo dhi
dlo dhi = change in box boundaries at end of run (distance units)
*scale* values = factor
factor = multiplicative factor for change in box length at end of run
*vel* value = V
V = change box length at this velocity (distance/time units),
effectively an engineering strain rate
*erate* value = R
R = engineering strain rate (1/time units)
*trate* value = R
R = true strain rate (1/time units)
*volume* value = none = adjust this dim to preserve volume of system
*wiggle* values = A Tp
A = amplitude of oscillation (distance units)
Tp = period of oscillation (time units)
*variable* values = v_name1 v_name2
v_name1 = variable with name1 for box length change as function of time
v_name2 = variable with name2 for change rate as function of time
*pressure* values = target gain
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
*pressure/mean* values = target gain
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
*xy*, *xz*, *yz* args = style value
style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure* or *erate/rescale*
*final* value = tilt
tilt = tilt factor at end of run (distance units)
*delta* value = dtilt
dtilt = change in tilt factor at end of run (distance units)
*vel* value = V
V = change tilt factor at this velocity (distance/time units),
effectively an engineering shear strain rate
*erate* value = R
R = engineering shear strain rate (1/time units)
*erate/rescale* value = R
R = engineering shear strain rate (1/time units)
*trate* value = R
R = true shear strain rate (1/time units)
*wiggle* values = A Tp
A = amplitude of oscillation (distance units)
Tp = period of oscillation (time units)
*variable* values = v_name1 v_name2
v_name1 = variable with name1 for tilt change as function of time
v_name2 = variable with name2 for change rate as function of time
*pressure* values = target gain
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
*erate/rescale* value = R
R = engineering shear strain rate (1/time units)
NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
*box* = style value
style = *volume* or *pressure*
@ -49,6 +86,15 @@ Syntax
.. parsed-literal::
*remap* value = *x* or *v* or *none*
x = remap coords of atoms in group into deforming box
v = remap velocities of atoms in group when they cross periodic boundaries
none = no remapping of x or v
*flip* value = *yes* or *no*
allow or disallow box flips when it becomes highly skewed
*units* value = *lattice* or *box*
lattice = distances are defined in lattice units
box = distances are defined in simulation box units
*couple* value = *none* or *xyz* or *xy* or *yz* or *xz*
couple pressure values of various dimensions
*vol/balance/p* value = *yes* or *no*
@ -57,7 +103,6 @@ Syntax
rate = maximum strain rate for pressure control
*normalize/pressure* value = *yes* or *no*
Modifies pressure controls such that the deviation in pressure is normalized by the target pressure
NOTE: All other keywords are documented by the :doc:`fix deform <fix_deform>` command
Examples
""""""""
@ -79,10 +124,26 @@ pressure-based controls implemented by this command.
All arguments described on the :doc:`fix deform <fix_deform>` doc page
also apply to this fix unless otherwise noted below. The rest of this
doc page explains the arguments specific to this fix. Note that a
page explains the arguments specific to this fix only. Note that a
simulation can define only a single deformation command: fix deform or
fix deform/pressure.
.. admonition:: Inconsistent trajectories due to image flags
:class: warning
When running long simulations while shearing the box or using a high
shearing rate, it is possible that the image flags used for storing
unwrapped atom positions will "wrap around". When LAMMPS is compiled
with the default settings, case image flags are limited to a range of
:math:`-512 \le i \le 511`, which will overflow when atoms starting
at zero image flag value have passed through a periodic box dimension
more than 512 times.
Changing the :ref:`size of LAMMPS integer types <size>` to the
"bigbig" setting can make this overflow much less likely, since it
increases the image flag value range to :math:`- 1,048,576 \le i \le
1\,048\,575`
----------
For the *x*, *y*, and *z* parameters, this is the meaning of the

22
doc/src/fix_electrode.rst
View File

@ -38,7 +38,7 @@ Syntax
*electrode/thermo* args = potential eta *temp* values
potential = electrode potential
charge = electrode charge
eta = reciprocal width of electrode charge smearing
eta = reciprocal width of electrode charge smearing (can be NULL if eta keyword is used)
*temp* values = T_v tau_v rng_v
T_v = temperature of thermo-potentiostat
tau_v = time constant of thermo-potentiostat
@ -287,8 +287,18 @@ The *fix_modify tf* option enables the Thomas-Fermi metallicity model
fix_modify ID tf type length voronoi
If this option is used parameters must be set for all atom types of the
electrode.
If this option is used, these two parameters must be set for
all atom types of the electrode:
* `tf` is the Thomas-Fermi length :math:`l_{TF}`
* `voronoi` is the Voronoi volume per atom in units of length cubed
Different types may have different `tf` and `voronoi` values.
The following self-energy term is then added for all electrode atoms:
.. math::
A_{ii} += \frac{1}{4 \pi \epsilon_0} \times \frac{4 \pi l_{TF}^2}{\mathrm{Voronoi volume}}
The *fix_modify timer* option turns on (off) additional timer outputs in the log
file, for code developers to track optimization.
@ -321,9 +331,11 @@ The global array has *N* rows and *2N+1* columns, where the fix manages
array, the elements are:
* array[I][1] = total charge that group *I* would have had *if it were
at 0 V applied potential* * array[I][2 to *N* + 1] = the *N* entries
at 0 V applied potential*
* array[I][2 to *N* + 1] = the *N* entries
of the *I*-th row of the electrode capacitance matrix (definition
follows) * array[I][*N* + 2 to *2N* + 1] = the *N* entries of the
follows)
* array[I][*N* + 2 to *2N* + 1] = the *N* entries of the
*I*-th row of the electrode elastance matrix (the inverse of the
electrode capacitance matrix)

10
doc/src/fix_heat_flow.rst
View File

@ -1,7 +1,7 @@
.. index:: fix heat/flow
fix heat/flow command
==========================
=====================
Syntax
""""""
@ -56,13 +56,19 @@ not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This pair style is part of the GRANULAR package. It is
only enabled if LAMMPS was built with that package.
See the :doc:`Build package <Build_package>` page for more info.
This fix requires that atoms store temperature and heat flow
as defined by the :doc:`fix property/atom <fix_property_atom>` command.
Related commands
""""""""""""""""
:doc:`pair granular <pair_granular>`, :doc:`fix property/atom <fix_property_atom>`
:doc:`pair granular <pair_granular>`,
:doc:`fix add/heat <fix_add_heat>`,
:doc:`fix property/atom <fix_property_atom>`
Default
"""""""

39
doc/src/fix_indent.rst
View File

@ -68,10 +68,10 @@ material or as an obstacle in a flow. Alternatively, it can be used as a
constraining wall around a simulation; see the discussion of the
*side* keyword below.
The *gstyle* geometry of the indenter can either be a sphere, a
cylinder, a cone, or a plane.
The *gstyle* keyword selects the geometry of the indenter and it can
either have the value of *sphere*, *cylinder*, *cone*, or *plane*\ .
A spherical indenter exerts a force of magnitude
A spherical indenter (*gstyle* = *sphere*) exerts a force of magnitude
.. math::
@ -82,13 +82,16 @@ distance from the atom to the center of the indenter, and *R* is the
radius of the indenter. The force is repulsive and F(r) = 0 for *r* >
*R*\ .
A cylindrical indenter exerts the same force, except that *r* is the
distance from the atom to the center axis of the cylinder. The
cylinder extends infinitely along its axis.
A cylindrical indenter (*gstyle* = *cylinder*) follows the same formula
for the force as a sphere, except that *r* is defined the distance
from the atom to the center axis of the cylinder. The cylinder extends
infinitely along its axis.
A conical indenter is similar to a cylindrical indenter except that it
has a finite length (between *lo* and *hi*), and that two different
radii (one at each end, *radlo* and *radhi*) can be defined.
.. versionadded:: 17April2024
A conical indenter (*gstyle* = *cone*) is similar to a cylindrical indenter
except that it has a finite length (between *lo* and *hi*), and that two
different radii (one at each end, *radlo* and *radhi*) can be defined.
Spherical, cylindrical, and conical indenters account for periodic
boundaries in two ways. First, the center point of a spherical
@ -101,15 +104,15 @@ or axis accounts for periodic boundaries. Both of these mean that an
indenter can effectively move through and straddle one or more
periodic boundaries.
A planar indenter is really an axis-aligned infinite-extent wall
exerting the same force on atoms in the system, where *R* is the
position of the plane and *r-R* is the distance from the plane. If
the *side* parameter of the plane is specified as *lo* then it will
indent from the lo end of the simulation box, meaning that atoms with
a coordinate less than the plane's current position will be pushed
towards the hi end of the box and atoms with a coordinate higher than
the plane's current position will feel no force. Vice versa if *side*
is specified as *hi*\ .
A planar indenter (*gstyle* = *plane*) behaves like an axis-aligned
infinite-extent wall with the same force expression on atoms in the
system as before, but where *R* is the position of the plane and *r-R*
is the distance of an from the plane. If the *side* parameter of the
plane is specified as *lo* then it will indent from the lo end of the
simulation box, meaning that atoms with a coordinate less than the
plane's current position will be pushed towards the hi end of the box
and atoms with a coordinate higher than the plane's current position
will feel no force. Vice versa if *side* is specified as *hi*\ .
Any of the 4 quantities defining a spherical indenter's geometry can
be specified as an equal-style :doc:`variable <variable>`, namely *x*,

180
doc/src/fix_rheo.rst Normal file
View File

@ -0,0 +1,180 @@
.. index:: fix rheo
fix rheo command
================
Syntax
""""""
.. parsed-literal::
fix ID group-ID rheo cut kstyle zmin keyword values...
* ID, group-ID are documented in :doc:`fix <fix>` command
* rheo = style name of this fix command
* cut = cutoff for the kernel (distance)
* kstyle = *quintic* or *RK0* or *RK1* or *RK2*
* zmin = minimal number of neighbors for reproducing kernels
* zero or more keyword/value pairs may be appended to args
* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or
*shift* or *rho/sum* or *density* or *self/mass* or *speed/sound*
.. parsed-literal::
*thermal* values = none, turns on thermal evolution
*interface/reconstruct* values = none, reconstructs interfaces with solid particles
*surface/detection* values = *sdstyle* *limit* *limit/splash*
*sdstyle* = *coordination* or *divergence*
*limit* = threshold for surface particles
*limit/splash* = threshold for splash particles
*shift* values = none, turns on velocity shifting
*rho/sum* values = none, uses the kernel to compute the density of particles
*self/mass* values = none, a particle uses its own mass in a rho summation
*density* values = *rho01*, ... *rho0N* (density)
*speed/sound* values = *cs0*, ... *csN* (velocity)
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all rheo 3.0 quintic 0 thermal density 0.1 0.1 speed/sound 10.0 1.0
fix 1 all rheo 3.0 RK1 10 shift surface/detection coordination 40
Description
"""""""""""
.. versionadded:: TBD
Perform time integration for RHEO particles, updating positions, velocities,
and densities. For a detailed breakdown of the integration timestep and
numerical details, see :ref:`(Palermo) <rheo_palermo>`. For an
overview of other features available in the RHEO package, see
:doc:`the RHEO howto <Howto_rheo>`.
The type of kernel is specified using *kstyle* and the cutoff is *cut*. Four
kernels are currently available. The *quintic* kernel is a standard quintic
spline function commonly used in SPH. The other options, *RK0*, *RK1*, and
*RK2*, are zeroth, first, and second order reproducing. To generate a
reproducing kernel, a particle must have sufficient neighbors inside the
kernel cutoff distance (a coordination number) to accurately calculate
moments. This threshold is set by *zmin*. If reproducing kernels are
requested but a particle has fewer neighbors, then it will revert to a
non-reproducing quintic kernel until it gains more neighbors.
To model temperature evolution, one must specify the *thermal* keyword,
define a separate instance of :doc:`fix rheo/thermal <fix_rheo_thermal>`,
and use atom style rheo/thermal.
By default, the density of solid RHEO particles does not evolve and forces
with fluid particles are calculated using the current velocity of the solid
particle. If the *interface/reconstruct* keyword is used, then the density
and velocity of solid particles are alternatively reconstructed for every
fluid-solid interaction to ensure no-slip and pressure-balanced boundaries.
This is done by estimating the location of the fluid-solid interface and
extrapolating fluid particle properties across the interface to calculate a
temporary apparent density and velocity for a solid particle. The numerical
details are the same as those described in
:ref:`(Palermo) <howto_rheo_palermo>` except there is an additional
restriction that the reconstructed solid density cannot be less than the
equilibrium density. This prevents fluid particles from sticking to solid
surfaces.
A modified form of Fickian particle shifting can be enabled with the
*shift* keyword. This effectively shifts particle positions to generate a
more uniform spatial distribution. Shifting currently does not consider the
type of a particle and therefore may be inappropriate in systems consisting
of multiple fluid phases.
In systems with free surfaces, the *surface/detection* keyword can be used
to classify the location of particles as being within the bulk fluid, on a
free surface, or isolated from other particles in a splash or droplet.
Shifting is then disabled in the normal direction away from the free surface
to prevent particles from diffusing away. Surface detection can also be used
to control surface-nucleated effects like oxidation when used in combination
with :doc:`fix rheo/oxidation <fix_rheo_oxidation>`. Surface detection is not
performed on solid bodies.
The *surface/detection* keyword takes three arguments: *sdstyle*, *limit*,
and *limit/splash*. The first, *sdstyle*, specifies whether surface particles
are identified using a coordination number (*coordination*) or the divergence
of the local particle positions (*divergence*). The threshold value for a
surface particle for either of these criteria is set by the numerical value
of *limit*. Additionally, if a particle's coordination number is too low,
i.e. if it has separated off from the bulk in a droplet, it is not possible
to define surfaces and the particle is classified as a splash. The coordination
threshold for this classification is set by the numerical value of
*limit/splash*.
By default, RHEO integrates particles' densities using a mass diffusion
equation. Alternatively, one can update densities every timestep by performing
a kernel summation of the masses of neighboring particles by specifying the *rho/sum*
keyword.
The *self/mass* keyword modifies the behavior of the density summation in *rho/sum*.
Typically, the density :math:`\rho` of a particle is calculated as the sum over neighbors
.. math::
\rho_i = \sum_{j} W_{ij} M_j
where :math:`W_{ij}` is the kernel, and :math:`M_j` is the mass of particle :math:`j`.
The *self/mass* keyword augments this expression by replacing :math:`M_j` with
:math:`M_i`. This may be useful in simulations of multiple fluid phases with large
differences in density, :ref:`(Hu) <fix_rheo_hu>`.
The *density* keyword is used to specify the equilibrium density of each of the N
particle types. It must be followed by N numerical values specifying each type's
equilibrium density *rho0*.
The *speed/sound* keyword is used to specify the speed of sound of each of the
N particle types. It must be followed by N numerical values specifying each type's
speed of sound *cs*.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix must be used with atom style rheo or rheo/thermal. This fix must
be used in conjunction with :doc:`fix rheo/pressure <fix_rheo_pressure>`.
and :doc:`fix rheo/viscosity <fix_rheo_viscosity>`. If the *thermal* setting
is used, there must also be an instance of
:doc:`fix rheo/thermal <fix_rheo_thermal>`. The fix group must be set to all.
Only one instance of fix rheo may be defined and it must be defined prior
to all other RHEO fixes in the input script.
This fix is part of the RHEO package. It is only enabled if LAMMPS was built
with that package. See the :doc:`Build package <Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo/viscosity <fix_rheo_viscosity>`,
:doc:`fix rheo/pressure <fix_rheo_pressure>`,
:doc:`fix rheo/thermal <fix_rheo_thermal>`,
:doc:`pair rheo <pair_rheo>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default
"""""""
*rho0* and *cs* are set to 1.0 for all atom types.
----------
.. _rheo_palermo:
**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, in preparation.
.. _fix_rheo_hu:
**(Hu)** Hu, and Adams J. Comp. Physics, 213, 844-861 (2006).

85
doc/src/fix_rheo_oxidation.rst Normal file
View File

@ -0,0 +1,85 @@
.. index:: fix rheo/oxidation
fix rheo/oxidation command
==========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID rheo/oxidation cut btype rsurf
* ID, group-ID are documented in :doc:`fix <fix>` command
* rheo/oxidation = style name of this fix command
* cut = maximum bond length (distance units)
* btype = type of bonds created
* rsurf = distance from surface to create bonds (distance units)
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all rheo/oxidation 1.5 2 0.0
fix 1 all rheo/oxidation 1.0 1 2.0
Description
"""""""""""
.. versionadded:: TBD
This fix dynamically creates bonds on the surface of fluids to
represent physical processes such as oxidation. It is intended
for use with bond style :doc:`bond rheo/shell <bond_rheo_shell>`.
Every timestep, particles check neighbors within a distance of *cut*.
This distance must be smaller than the kernel length defined in
:doc:`fix rheo <fix_rheo>`. Bonds of type *btype* are created between
a fluid particle and either a fluid or solid neighbor. The fluid particles
must also be on the fluid surface, or within a distance of *rsurf* from
the surface. This process is further described in
:ref:`(Clemmer) <howto_rheo_clemmer2>`.
If used in conjunction with solid bodies, such as those generated
by the *react* option of :doc:`fix rheo/thermal <fix_rheo_thermal>`,
it is recommended to use a :doc:`hybrid bond style <bond_hybrid>`
with different bond types for solid and oxide bonds.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix must be used with the bond style :doc:`rheo/shell <bond_rheo_shell>`
and :doc:`fix rheo <fix_rheo>` with surface detection enabled.
This fix is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>`
page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`bond rheo/shell <bond_rheo_shell>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default
"""""""
none
----------
.. _howto_rheo_clemmer2:
**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

106
doc/src/fix_rheo_pressure.rst Normal file
View File

@ -0,0 +1,106 @@
.. index:: fix rheo/pressure
fix rheo/pressure command
=========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID rheo/pressure type1 pstyle1 args1 ... typeN pstyleN argsN
* ID, group-ID are documented in :doc:`fix <fix>` command
* rheo/pressure = style name of this fix command
* one or more types and pressure styles must be appended
* types = lists of types (see below)
* pstyle = *linear* or *taitwater* or *cubic*
.. parsed-literal::
*linear* args = none
*taitwater* args = none
*cubic* args = cubic prefactor :math:`A_3` (pressure/density\^2)
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all rheo/pressure * linear
fix 1 all rheo/pressure 1 linear 2 cubic 10.0
Description
"""""""""""
.. versionadded:: TBD
This fix defines a pressure equation of state for RHEO particles. One can
define different equations of state for different atom types. An equation
must be specified for every atom type.
One first defines the atom *types*. A wild-card asterisk can be used in place
of or in conjunction with the *types* argument to set the coefficients for
multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*"
or "m\*n". If :math:`N` is the number of atom types, then an asterisk with
no numeric values means all types from 1 to :math:`N`. A leading asterisk
means all types from 1 to n (inclusive). A trailing asterisk means all types
from m to :math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
The *types* definition is followed by the pressure style, *pstyle*. Current
options *linear*, *taitwater*, and *cubic*. Style *linear* is a linear
equation of state with a particle pressure :math:`P` calculated as
.. math::
P = c (\rho - \rho_0)
where :math:`c` is the speed of sound, :math:`\rho_0` is the equilibrium density,
and :math:`\rho` is the current density of a particle. The numerical values of
:math:`c` and :math:`\rho_0` are set in :doc:`fix rheo <fix_rheo>`. Style *cubic*
is a cubic equation of state which has an extra argument :math:`A_3`,
.. math::
P = c ((\rho - \rho_0) + A_3 (\rho - \rho_0)^3) .
Style *taitwater* is Tait's equation of state:
.. math::
P = \frac{c^2 \rho_0}{7} \biggl[\left(\frac{\rho}{\rho_0}\right)^{7} - 1\biggr].
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix must be used with an atom style that includes density
such as atom_style rheo or rheo/thermal. This fix must be used in
conjunction with :doc:`fix rheo <fix_rheo>`. The fix group must be
set to all. Only one instance of fix rheo/pressure can be defined.
This fix is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>`
page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`pair rheo <pair_rheo>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default
"""""""
none

128
doc/src/fix_rheo_thermal.rst Normal file
View File

@ -0,0 +1,128 @@
.. index:: fix rheo/thermal
fix rheo/thermal command
========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID rheo/thermal attribute values ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* rheo/thermal = style name of this fix command
* one or more attributes may be appended
* attribute = *conductivity* or *specific/heat* or *latent/heat* or *Tfreeze* or *react*
.. parsed-literal::
*conductivity* args = types style args
types = lists of types (see below)
style = *constant*
*constant* arg = conductivity (power/temperature)
*specific/heat* args = types style args
types = lists of types (see below)
style = *constant*
*constant* arg = specific heat (energy/(mass*temperature))
*latent/heat* args = types style args
types = lists of types (see below)
style = *constant*
*constant* arg = latent heat (energy/mass)
*Tfreeze* args = types style args
types = lists of types (see below)
style = *constant*
*constant* arg = freezing temperature (temperature)
*react* args = cut type
cut = maximum bond distance
type = bond type
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all rheo/thermal conductivity * constant 1.0 specific/heat * constant 1.0 Tfreeze * constant 1.0
fix 1 all rheo/pressure conductivity 1*2 constant 1.0 conductivity 3*4 constant 2.0 specific/heat * constant 1.0
Description
"""""""""""
.. versionadded:: TBD
This fix performs time integration of temperature for atom style rheo/thermal.
In addition, it defines multiple thermal properties of particles and handles
melting/solidification, if applicable. For more details on phase transitions
in RHEO, see :doc:`the RHEO howto <Howto_rheo>`.
Note that the temperature of a particle is always derived from the energy.
This implies the *temperature* attribute of :doc:`the set command <set>` does
not affect particles. Instead, one should use the *sph/e* attribute.
For each atom type, one can define expressions for the *conductivity*,
*specific/heat*, *latent/heat*, and critical temperature (*Tfreeze*).
The conductivity and specific heat must be defined for all atom types.
The latent heat and critical temperature are optional. However, a
critical temperature must be defined to specify a latent heat.
Note, if shifting is turned on in :doc:`fix rheo <fix_rheo>`, the gradient
of the energy is used to shift energies. This may be inappropriate in systems
with multiple atom types with different specific heats.
For each property, one must first define a list of atom types. A wild-card
asterisk can be used in place of or in conjunction with the *types* argument
to set the coefficients for multiple pairs of atom types. This takes the
form "\*" or "\*n" or "m\*" or "m\*n". If :math:`N` is the number of atom
types, then an asterisk with no numeric values means all types from 1 to
:math:`N`. A leading asterisk means all types from 1 to n (inclusive).
A trailing asterisk means all types from m to :math:`N` (inclusive). A
middle asterisk means all types from m to n (inclusive).
The *types* definition for each property is followed by the style. Currently,
the only option is *constant*. Style *constant* simply applies a constant value
of respective property to each particle of the assigned type.
The *react* keyword controls whether bonds are created/deleted when particles
transition between a fluid and solid state. This option only applies to atom
types that have a defined value of *Tfreeze*. When a fluid particle's
temperature drops below *Tfreeze*, bonds of type *btype* are created between
nearby solid particles within a distance of *cut*. The particle's status also
swaps to a solid state. When a solid particle's temperature rises above
*Tfreeze*, all bonds of type *btype* are broken and the particle's status swaps
to a fluid state.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix must be used with an atom style that includes temperature,
heatflow, and conductivity such as atom_style rheo/thermal This fix
must be used in conjunction with :doc:`fix rheo <fix_rheo>` with the
*thermal* setting. The fix group must be set to all. Only one
instance of fix rheo/pressure can be defined.
This fix is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>`
page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`pair rheo <pair_rheo>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`,
:doc:`fix add/heat <fix_add_heat>`
Default
"""""""
none

117
doc/src/fix_rheo_viscosity.rst Normal file
View File

@ -0,0 +1,117 @@
.. index:: fix rheo/viscosity
fix rheo/viscosity command
==========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID rheo/viscosity type1 pstyle1 args1 ... typeN pstyleN argsN
* ID, group-ID are documented in :doc:`fix <fix>` command
* rheo/viscosity = style name of this fix command
* one or more types and viscosity styles must be appended
* types = lists of types (see below)
* vstyle = *constant* or *power*
.. parsed-literal::
*constant* args = *eta*
*eta* = viscosity
*power* args = *eta*, *gd0*, *K*, *n*
*eta* = viscosity
*gd0* = critical strain rate
*K* = consistency index
*n* = power-law exponent
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all rheo/viscosity * constant 1.0
fix 1 all rheo/viscosity 1 constant 1.0 2 power 0.1 5e-4 0.001 0.5
Description
"""""""""""
.. versionadded:: TBD
This fix defines a viscosity for RHEO particles. One can define different
viscosities for different atom types, but a viscosity must be specified for
every atom type.
One first defines the atom *types*. A wild-card asterisk can be used in place
of or in conjunction with the *types* argument to set the coefficients for
multiple pairs of atom types. This takes the form "\*" or "\*n" or "m\*"
or "m\*n". If :math:`N` is the number of atom types, then an asterisk with
no numeric values means all types from 1 to :math:`N`. A leading asterisk
means all types from 1 to n (inclusive). A trailing asterisk means all types
from m to :math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
The *types* definition is followed by the viscosity style, *vstyle*. Two
options are available, *constant* and *power*. Style *constant* simply
applies a constant value of the viscosity *eta* to each particle of the
assigned type. Style *power* is a Hershchel-Bulkley constitutive equation
for the stress :math:`\tau`
.. math::
\tau = \left(\frac{\tau_0}{\dot{\gamma}} + K \dot{\gamma}^{n - 1}\right) \dot{\gamma}, \tau \ge \tau_0
where :math:`\dot{\gamma}` is the strain rate and :math:`\tau_0` is the critical
yield stress, below which :math:`\dot{\gamma} = 0.0`. To avoid divergences, this
expression is regularized by defining a critical strain rate *gd0*. If the local
strain rate on a particle falls below this limit, a constant viscosity of *eta*
is assigned. This implies a value of
.. math::
\tau_0 = \eta \dot{\gamma}_0 - K \dot{\gamma}_0^N
as further discussed in :ref:`(Palermo) <rheo_palermo2>`.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`.
None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix. No global or per-atom quantities are stored
by this fix for access by various :doc:`output commands <Howto_output>`.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
This fix must be used with an atom style that includes viscosity
such as atom_style rheo or rheo/thermal. This fix must be used in
conjunction with :doc:`fix rheo <fix_rheo>`. The fix group must be
set to all. Only one instance of fix rheo/viscosity can be defined.
This fix is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`pair rheo <pair_rheo>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default
"""""""
none
----------
.. _rheo_palermo2:
**(Palermo)** Palermo, Wolf, Clemmer, O'Connor, in preparation.

26
doc/src/fix_store_force.rst
View File

@ -23,11 +23,12 @@ Examples
Description
"""""""""""
Store the forces on atoms in the group at the point during each
timestep when the fix is invoked, as described below. This is useful
for storing forces before constraints or other boundary conditions are
computed which modify the forces, so that unmodified forces can be
:doc:`written to a dump file <dump>` or accessed by other :doc:`output commands <Howto_output>` that use per-atom quantities.
Store the forces on atoms in the group at the point during each timestep
when the fix is invoked, as described below. This is useful for storing
forces before constraints or other boundary conditions are computed
which modify the forces, so that unmodified forces can be :doc:`written
to a dump file <dump>` or accessed by other :doc:`output commands
<Howto_output>` that use per-atom quantities.
This fix is invoked at the point in the velocity-Verlet timestepping
immediately after :doc:`pair <pair_style>`, :doc:`bond <bond_style>`,
@ -36,12 +37,13 @@ immediately after :doc:`pair <pair_style>`, :doc:`bond <bond_style>`,
forces have been calculated. It is the point in the timestep when
various fixes that compute constraint forces are calculated and
potentially modify the force on each atom. Examples of such fixes are
:doc:`fix shake <fix_shake>`, :doc:`fix wall <fix_wall>`, and :doc:`fix indent <fix_indent>`.
:doc:`fix shake <fix_shake>`, :doc:`fix wall <fix_wall>`, and :doc:`fix
indent <fix_indent>`.
.. note::
The order in which various fixes are applied which operate at
the same point during the timestep, is the same as the order they are
The order in which various fixes are applied which operate at the
same point during the timestep, is the same as the order they are
specified in the input script. Thus normally, if you want to store
per-atom forces due to force field interactions, before constraints
are applied, you should list this fix first within that set of fixes,
@ -52,8 +54,9 @@ potentially modify the force on each atom. Examples of such fixes are
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`. None of the :doc:`fix_modify <fix_modify>` options
are relevant to this fix.
No information about this fix is written to :doc:`binary restart files
<restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix.
This fix produces a per-atom array which can be accessed by various
:doc:`output commands <Howto_output>`. The number of columns for each
@ -61,7 +64,8 @@ atom is 3, and the columns store the x,y,z forces on each atom. The
per-atom values be accessed on any timestep.
No parameter of this fix can be used with the *start/stop* keywords of
the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minimization <minimize>`.
the :doc:`run <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""

27
doc/src/pair_coul.rst
View File

@ -2,6 +2,8 @@
.. index:: pair_style coul/cut/gpu
.. index:: pair_style coul/cut/kk
.. index:: pair_style coul/cut/omp
.. index:: pair_style coul/cut/global
.. index:: pair_style coul/cut/global/omp
.. index:: pair_style coul/debye
.. index:: pair_style coul/debye/gpu
.. index:: pair_style coul/debye/kk
@ -11,8 +13,6 @@
.. index:: pair_style coul/dsf/kk
.. index:: pair_style coul/dsf/omp
.. index:: pair_style coul/exclude
.. index:: pair_style coul/cut/global
.. index:: pair_style coul/cut/global/omp
.. index:: pair_style coul/long
.. index:: pair_style coul/long/omp
.. index:: pair_style coul/long/kk
@ -33,6 +33,11 @@ pair_style coul/cut command
Accelerator Variants: *coul/cut/gpu*, *coul/cut/kk*, *coul/cut/omp*
pair_style coul/cut/global command
==================================
Accelerator Variants: *coul/cut/omp*
pair_style coul/debye command
=============================
@ -46,11 +51,6 @@ Accelerator Variants: *coul/dsf/gpu*, *coul/dsf/kk*, *coul/dsf/omp*
pair_style coul/exclude command
===============================
pair_style coul/cut/global command
==================================
Accelerator Variants: *coul/cut/omp*
pair_style coul/long command
============================
@ -79,16 +79,17 @@ pair_style tip4p/long command
Accelerator Variants: *tip4p/long/omp*
Syntax
""""""
.. code-block:: LAMMPS
pair_style coul/cut cutoff
pair_style coul/cut/global cutoff
pair_style coul/debye kappa cutoff
pair_style coul/dsf alpha cutoff
pair_style coul/exclude cutoff
pair_style coul/cut/global cutoff
pair_style coul/long cutoff
pair_style coul/wolf alpha cutoff
pair_style coul/streitz cutoff keyword alpha
@ -152,6 +153,11 @@ the 2 atoms, and :math:`\epsilon` is the dielectric constant which can be set
by the :doc:`dielectric <dielectric>` command. The cutoff :math:`r_c` truncates
the interaction distance.
Pair style *coul/cut/global* computes the same Coulombic interactions
as style *coul/cut* except that it allows only a single global cutoff
and thus makes it compatible for use in combination with long-range
coulomb styles in :doc:`hybrid pair styles <pair_hybrid>`.
----------
Style *coul/debye* adds an additional exp() damping factor to the
@ -262,11 +268,6 @@ Streitz-Mintmire parameterization for the material being modeled.
----------
Pair style *coul/cut/global* computes the same Coulombic interactions
as style *coul/cut* except that it allows only a single global cutoff
and thus makes it compatible for use in combination with long-range
coulomb styles in :doc:`hybrid pair styles <pair_hybrid>`.
Pair style *coul/exclude* computes Coulombic interactions like *coul/cut*
but **only** applies them to excluded pairs using a scaling factor
of :math:`\gamma - 1.0` with :math:`\gamma` being the factor assigned

102
doc/src/pair_rheo.rst Normal file
View File

@ -0,0 +1,102 @@
.. index:: pair_style rheo
pair_style rheo command
=======================
Syntax
""""""
.. code-block:: LAMMPS
pair_style rheo cutoff keyword values
* cutoff = global cutoff for kernel (distance units)
* zero or more keyword/value pairs may be appended to args
* keyword = *rho/damp* or *artificial/visc* or *harmonic/means*
.. parsed-literal::
*rho/damp* args = density damping prefactor :math:`\xi`
*artificial/visc* args = artificial viscosity prefactor :math:`\zeta`
*harmonic/means* args = none
Examples
""""""""
.. code-block:: LAMMPS
pair_style rheo 3.0 rho/damp 1.0 artificial/visc 2.0
pair_coeff * *
Description
"""""""""""
.. versionadded:: TBD
Pair style *rheo* computes pressure and viscous forces between particles
in the :doc:`rheo package <Howto_rheo>`. If thermal evolution is turned
on in :doc:`fix rheo <fix_rheo>`, then the pair style also calculates
heat exchanged between particles.
The *artificial/viscosity* keyword is used to specify the magnitude
:math:`\zeta` of an optional artificial viscosity contribution to forces.
This factor can help stabilize simulations by smoothing out small length
scale variations in velocity fields. Artificial viscous forces typically
are only exchanged by fluid particles. However, if interfaces are not
reconstructed in fix rheo, fluid particles will also exchange artificial
viscous forces with solid particles to improve stability.
The *rho/damp* keyword is used to specify the magnitude :math:`\xi` of
an optional pairwise damping term between the density of particles. This
factor can help stabilize simulations by smoothing out small length
scale variations in density fields. However, in systems that develop
a density gradient in equilibrium (e.g. in a hydrostatic column underlying
gravity), this option may be inappropriate.
If particles have different viscosities or conductivities, the
*harmonic/means* keyword changes how they are averaged before calculating
pairwise forces or heat exchanges. By default, an arithmetic averaged is
used, however, a harmonic mean may improve stability in systems with multiple
fluid phases with large disparities in viscosities.
No coefficients are defined for each pair of atoms types via the
:doc:`pair_coeff <pair_coeff>` command as in the examples
above.
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This style does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options.
This style does not write information to :doc:`binary restart files <restart>`.
Thus, you need to re-specify the pair_style and pair_coeff commands in an input
script that reads a restart file.
This style can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the *inner*,
*middle*, *outer* keywords.
Restrictions
""""""""""""
This fix is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`fix rheo/pressure <fix_rheo_pressure>`,
:doc:`fix rheo/thermal <fix_rheo_thermal>`,
:doc:`fix rheo/viscosity <fix_rheo_viscosity>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default
"""""""
Density damping and artificial viscous forces are not calculated.
Arithmetic means are used for mixing particle properties.

112
doc/src/pair_rheo_solid.rst Normal file
View File

@ -0,0 +1,112 @@
.. index:: pair_style rheo/solid
pair_style rheo/solid command
=============================
Syntax
""""""
.. code-block:: LAMMPS
pair_style rheo/solid
Examples
""""""""
.. code-block:: LAMMPS
pair_style rheo/solid
pair_coeff * * 1.0 1.5 1.0
Description
"""""""""""
.. versionadded:: TBD
Style *rheo/solid* is effectively a copy of pair style
:doc:`bpm/spring <pair_bpm_spring>` except it only applies forces
between solid RHEO particles, determined by checking the status of
each pair of neighboring particles before calculating forces.
The style computes pairwise forces with the formula
.. math::
F = k (r - r_c)
where :math:`k` is a stiffness and :math:`r_c` is the cutoff length.
An additional damping force is also applied to interacting
particles. The force is proportional to the difference in the
normal velocity of particles
.. math::
F_D = - \gamma w (\hat{r} \bullet \vec{v})
where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the
displacement normal vector, :math:`\vec{v}` is the velocity difference
between the two particles, and :math:`w` is a smoothing factor.
This smoothing factor is constructed such that damping forces go to zero
as particles come out of contact to avoid discontinuities. It is
given by
.. math::
w = 1.0 - \left( \frac{r}{r_c} \right)^8 .
The following coefficients must be defined for each pair of atom 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, or by mixing as described below:
* :math:`k` (force/distance units)
* :math:`r_c` (distance units)
* :math:`\gamma` (force/velocity units)
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
For atom type pairs I,J and I != J, the A coefficient and cutoff
distance for this pair style can be mixed. A is always mixed via a
*geometric* rule. The cutoff is mixed according to the pair_modify
mix value. The default mix value is *geometric*\ . See the
"pair_modify" command for details.
This pair style does not support the :doc:`pair_modify <pair_modify>`
shift option, since the pair interaction goes to 0.0 at the cutoff.
The :doc:`pair_modify <pair_modify>` table and tail options are not
relevant for this pair style.
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 can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the
*inner*, *middle*, *outer* keywords.
----------
Restrictions
""""""""""""
This pair style is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`fix rheo/thermal <fix_rheo_thermal>`,
:doc:`pair bpm/spring <pair_bpm_spring>`
Default
"""""""
none

2
doc/src/pair_style.rst
View File

@ -337,6 +337,8 @@ accelerated styles exist.
* :doc:`reaxff <pair_reaxff>` - ReaxFF potential
* :doc:`rebo <pair_airebo>` - Second generation REBO potential of Brenner
* :doc:`rebomos <pair_rebomos>` - REBOMoS potential for MoS2
* :doc:`rheo <pair_rheo>` - fluid interactions in RHEO package
* :doc:`rheo/solid <pair_rheo_solid>` - solid interactions in RHEO package
* :doc:`resquared <pair_resquared>` - Everaers RE-Squared ellipsoidal potential
* :doc:`saip/metal <pair_saip_metal>` - Interlayer potential for hetero-junctions formed with hexagonal 2D materials and metal surfaces
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>` - Smoothed dissipative particle dynamics for water at isothermal conditions

6
doc/src/read_data.rst
View File

@ -12,7 +12,7 @@ Syntax
* file = name of data file to read in
* zero or more keyword/arg pairs may be appended
* keyword = *add* or *offset* or *shift* or *extra/atom/types* or *extra/bond/types* or *extra/angle/types* or *extra/dihedral/types* or *extra/improper/types* or *extra/bond/per/atom* or *extra/angle/per/atom* or *extra/dihedral/per/atom* or *extra/improper/per/atom* or *group* or *nocoeff* or *fix*
* keyword = *add* or *offset* or *shift* or *extra/atom/types* or *extra/bond/types* or *extra/angle/types* or *extra/dihedral/types* or *extra/improper/types* or *extra/bond/per/atom* or *extra/angle/per/atom* or *extra/dihedral/per/atom* or *extra/improper/per/atom* or *extra/special/per/atom* or *group* or *nocoeff* or *fix*
.. parsed-literal::
@ -859,6 +859,10 @@ of analysis.
- atom-ID molecule-ID atom-type x y z
* - peri
- atom-ID atom-type volume density x y z
* - rheo
- atom-ID atom-type status rho x y z
* - rheo/thermal
- atom-ID atom-type status rho energy x y z
* - smd
- atom-ID atom-type molecule volume mass kradius cradius x0 y0 z0 x y z
* - sph

22
doc/src/reset_atoms.rst
View File

@ -32,7 +32,7 @@ Syntax
.. code-block:: LAMMPS
reset atoms mol group-ID keyword value ...
reset_atoms mol group-ID keyword value ...
* group-ID = ID of group of atoms whose molecule IDs will be reset
* zero or more keyword/value pairs can be appended
@ -66,16 +66,16 @@ Description
.. versionadded:: 22Dec2022
The *reset_atoms* command resets the values of a specified atom
property. In contrast to the set command, it does this in a
property. In contrast to the *set* command, it does this in a
collective manner which resets the values for many atoms in a
self-consistent way. This is often useful when the simulated system
has undergone significant modifications like adding or removing atoms
or molecules, joining data files, changing bonds, or large-scale
self-consistent way. This command is often useful when the simulated
system has undergone significant modifications like adding or removing
atoms or molecules, joining data files, changing bonds, or large-scale
diffusion.
The new values can be thought of as a *reset*, similar to values atoms
would have if a new data file were being read or a new simulation
performed. Note that the set command also resets atom properties to
performed. Note that the *set* command also resets atom properties to
new values, but it treats each atom independently.
The *property* setting can be *id* or *image* or *mol*. For *id*, the
@ -90,7 +90,7 @@ keyword/value settings are given below.
----------
*Property id*
Property: *id*
Reset atom IDs for the entire system, including all the global IDs
stored for bond, angle, dihedral, improper topology data. This will
@ -146,7 +146,7 @@ processor have consecutive IDs, as the :doc:`create_atoms
----------
*Property image*
Property: *image*
Reset the image flags of atoms so that at least one atom in each
molecule has an image flag of 0. Molecular topology is respected so
@ -191,7 +191,7 @@ flags.
----------
*Property mol*
Property: *mol*
Reset molecule IDs for a specified group of atoms based on current
bond connectivity. This will typically create a new set of molecule
@ -203,7 +203,7 @@ For purposes of this operation, molecules are identified by the current
bond connectivity in the system, which may or may not be consistent with
the current molecule IDs. A molecule in this context is a set of atoms
connected to each other with explicit bonds. The specific algorithm
used is the one of :doc:`compute fragment/atom <compute_cluster_atom>`
used is the one of :doc:`compute fragment/atom <compute_cluster_atom>`.
Once the molecules are identified and a new molecule ID computed for
each, this command will update the current molecule ID for all atoms in
the group with the new molecule ID. Note that if the group excludes
@ -266,7 +266,7 @@ The *image* property can only be used when the atom style supports bonds.
Related commands
""""""""""""""""
:doc:`compute fragment/atom <compute_cluster_atom>`
:doc:`compute fragment/atom <compute_cluster_atom>`,
:doc:`fix bond/react <fix_bond_react>`,
:doc:`fix bond/create <fix_bond_create>`,
:doc:`fix bond/break <fix_bond_break>`,

6
doc/src/set.rst
View File

@ -120,6 +120,8 @@ Syntax
*angle* value = numeric angle type or angle type label, for all angles between selected atoms
*dihedral* value = numeric dihedral type or dihedral type label, for all dihedrals between selected atoms
*improper* value = numeric improper type or improper type label, for all impropers between selected atoms
*rheo/rho* value = density of RHEO particles (mass/distance\^3)
*rheo/status* value = status or phase of RHEO particles (unitless)
*sph/e* value = energy of SPH particles (need units)
value can be an atom-style variable (see below)
*sph/cv* value = heat capacity of SPH particles (need units)
@ -506,6 +508,10 @@ by the *bond types* (\ *angle types*, etc) field in the header of the
data file read by the :doc:`read_data <read_data>` command. These
keywords do not allow use of an atom-style variable.
Keywords *rheo/rho* and *rheo/status* set the density and the status of
rheo particles. In particular, one can only set the phase in the status
as described by the :doc:`RHEO howto page <Howto_rheo>`.
Keywords *sph/e*, *sph/cv*, and *sph/rho* set the energy, heat capacity,
and density of smoothed particle hydrodynamics (SPH) particles. See
`this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in LAMMPS.

4
doc/utils/sphinx-config/conf.py.in
View File

@ -41,7 +41,7 @@ sys.path.append(os.path.join(LAMMPS_DOC_DIR, 'utils', 'sphinx-config', '_themes'
# -- General configuration ------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
needs_sphinx = '5.2.0'
needs_sphinx = '5.3.0'
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
@ -69,7 +69,7 @@ images_config = {
templates_path = ['_templates']
# The suffix of source filenames.
source_suffix = '.rst'
source_suffix = {'.rst': 'restructuredtext'}
# The encoding of source files.
#source_encoding = 'utf-8-sig'

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

@ -393,6 +393,7 @@ buf
builtin
Bulacu
Bulatov
Bulkley
Bureekaew
burlywood
Bussi
@ -564,6 +565,7 @@ cond
conda
Conda
Condens
conductivities
conf
config
configfile
@ -1440,6 +1442,7 @@ henrich
Henrich
Hermitian
Herrmann
Hershchel
Hertizian
hertzian
Hertzsch
@ -1831,6 +1834,7 @@ Kspace
KSpace
KSpaceStyle
Kspring
kstyle
kT
kTequil
kth
@ -2271,6 +2275,7 @@ modelled
modelling
Modelling
Modine
modularity
moduli
mofff
MOFFF
@ -2488,6 +2493,7 @@ Neumann
Nevent
nevery
Nevery
Nevins
newfile
Newns
newtype
@ -3067,6 +3073,7 @@ quatw
queryargs
Queteschiner
quickmin
quintic
qw
qx
qy
@ -3078,6 +3085,7 @@ radialscreenedspin
radialspin
radian
radians
radiative
radj
Rafferty
rahman
@ -3181,6 +3189,7 @@ rg
Rg
Rhaphson
Rhe
rheo
rheological
rheology
rhodo
@ -3268,6 +3277,7 @@ rsort
rsq
rst
rstyle
rsurf
Rubensson
Rubia
Rud
@ -3651,6 +3661,7 @@ Telsa
tempCorrCoeff
templated
Templeton
Tencer
Tequil
ters
tersoff
@ -3997,6 +4008,7 @@ Vries
Vsevolod
Vsmall
Vstream
vstyle
vtarget
vtk
VTK

1
examples/PACKAGES/electrode/madelung/.gitignore vendored
View File

@ -1,2 +1,3 @@
*.csv
*.txt
*.lammpstrj

18
examples/PACKAGES/electrode/madelung/eval.py
View File

@ -17,14 +17,22 @@ q_ref = float(ref_line[3])
inv11_ref = float(ref_line[4])
inv12_ref = float(ref_line[5])
b1_ref = float(ref_line[6])
felec1_ref = float(ref_line[8])
felyt1_ref = float(ref_line[10])
press_ref = float(ref_line[12])
# out.csv
with open(sys.argv[2]) as f:
out_line = f.readlines()[-1].split(", ")
e_out = float(out_line[0])
q_out = float(out_line[1])
press_out = float(out_line[2])
out_lines = [("energy", e_ref, e_out), ("charge", q_ref, q_out)]
out_lines = [
("energy", e_ref, e_out),
("charge", q_ref, q_out),
("pressure", press_ref, press_out),
]
# vec.csv
vec_file = "vec.csv"
@ -44,6 +52,14 @@ if op.isfile(inv_file):
inv12_out = float(inv_line[1])
out_lines.append(("inv11", inv11_ref, inv11_out))
# forces.lammpstrj
force_file = "forces.lammpstrj"
with open(force_file) as f:
lines = f.readlines()[9:]
for name, i, f_ref in [("felec1", "1", felec1_ref), ("felyt1", "3", felyt1_ref)]:
f_out = next(float(y[3]) for x in lines if (y := x.split()) and y[0] == i)
out_lines.append((name, f_ref, f_out))
lines = []
for label, ref, out in out_lines:
error = rel_error(out, ref)

2
examples/PACKAGES/electrode/madelung/in.eta_cg
View File

@ -8,7 +8,7 @@ thermo_style custom step pe c_qbot c_qtop
fix feta all property/atom d_eta ghost on
set group bot d_eta 0.5
set group top d_eta 3.0
fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta algo cg 1e-6
fix conp bot electrode/conp 0 NULL couple top 1 symm on eta d_eta algo cg 1e-6
run 0

2
examples/PACKAGES/electrode/madelung/in.eta_mix
View File

@ -8,7 +8,7 @@ thermo_style custom step pe c_qbot c_qtop
fix feta all property/atom d_eta ghost on
set group bot d_eta 0.5
set group top d_eta 3.0
fix conp bot electrode/conp 0 2 couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
fix conp bot electrode/conp 0 NULL couple top 1 symm on eta d_eta write_inv inv.csv write_vec vec.csv
run 0

98
examples/PACKAGES/electrode/madelung/plate_cap.py
View File

@ -1,12 +1,17 @@
#!/usr/bin/env python3
import time
import numpy as np
from scipy.special import erf
SQRT2 = np.sqrt(2)
SQRTPI_INV = 1 / np.sqrt(np.pi)
COULOMB = 332.06371 # Coulomb constant in Lammps 'real' units
QE2F = 23.060549
NKTV2P = 68568.415 # pressure in 'real' units
LENGTH = 10000 # convergence parameter
LZ = 20
def lattice(length):
@ -26,6 +31,25 @@ def b_element(r, q, eta):
return q * erf(eta * r) / r
def force_gauss(r, qq, eta):
etar = eta * r
return (qq / np.square(r)) * (
erf(etar) - 2 * etar * SQRTPI_INV * np.exp(-np.square(etar))
)
def force_point(r, qq):
return qq / np.square(r)
def force_component(dx, d, qq, eta=None):
if eta:
return np.sum(dx / d * force_gauss(d, qq, eta))
else:
return np.sum(dx / d * force_point(d, qq))
time_start = time.perf_counter()
a = 1 # nearest neighbor distance i.e. lattice constant / sqrt(2)
x_elec = [-2, 2]
x_elyt = [-1, 1]
@ -36,8 +60,20 @@ v = np.array([-0.5, 0.5]) * (QE2F / COULOMB)
# distances to images within electrode and to opposite electrode
distances = a * np.linalg.norm(lattice(LENGTH), axis=1)
opposite_distances = np.sqrt(np.square(distances) + distance_plates**2)
image_distances = []
for x in x_elec:
image_distances.append([])
for y in x_elyt:
image_distances[-1].append(np.sqrt(np.square(distances) + np.abs(y - x) ** 2))
image_elyt_distances = [[None for _ in range(len(x_elyt))] for _ in range(len(x_elyt))]
for i, (xi, qi) in enumerate(zip(x_elyt, q_elyt)):
for j, (xj, qj) in list(enumerate(zip(x_elyt, q_elyt)))[i + 1 :]:
image_elyt_distances[i][j] = np.sqrt(
np.square(distances) + np.abs(xj - xi) ** 2
)
for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
# for name, eta_elec in [("", [2.0, 2.0])]:
eta_mix = np.prod(eta_elec) / np.sqrt(np.sum(np.square(eta_elec)))
# self interaction and within original box
A_11 = np.sqrt(2 / np.pi) * eta_elec[0]
@ -55,22 +91,18 @@ for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
# electrode-electrolyte interaction
b = []
for x, eta in zip(x_elec, eta_elec):
for i, (x, eta) in enumerate(zip(x_elec, eta_elec)):
bi = 0
for y, q in zip(x_elyt, q_elyt):
d = abs(y - x)
bi += b_element(d, q, eta)
image_distances = np.sqrt(np.square(distances) + d**2)
bi += 4 * np.sum(b_element(image_distances, q, eta))
for j, (y, q) in enumerate(zip(x_elyt, q_elyt)):
bi += b_element(np.abs(y - x), q, eta)
bi += 4 * np.sum(b_element(image_distances[i][j], q, eta))
b.append(bi)
b = np.array(b)
# electrolyte-electrolyte energy
elyt_11 = 4 * np.sum(1 / distances)
distance_elyt = x_elyt[1] - x_elyt[0]
elyt_12 = 1 / distance_elyt + 4 * np.sum(
1 / np.sqrt(np.square(distances) + distance_elyt**2)
)
elyt_12 = 1 / distance_elyt + 4 * np.sum(1 / image_elyt_distances[0][1])
elyt = np.array([[elyt_11, elyt_12], [elyt_12, elyt_11]])
energy_elyt = 0.5 * np.dot(q_elyt, np.dot(elyt, q_elyt))
@ -78,9 +110,48 @@ for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
q = np.dot(inv, v - b)
energy = COULOMB * (0.5 * np.dot(q, np.dot(A, q)) + np.dot(b, q) + energy_elyt)
# forces in out-of-plane direction
f_elec = np.zeros(len(x_elec))
f_elyt = np.zeros(len(x_elyt))
# electrode-electrode
dx = x_elec[1] - x_elec[0]
fij_box = force_component(dx, np.abs(dx), q[0] * q[1], eta_mix)
fij_img = 4 * force_component(dx, opposite_distances, q[0] * q[1], eta_mix)
f_elec[0] -= fij_box + fij_img
f_elec[1] += fij_box + fij_img
# electrode-electrolyte
for i, (xi, qi, etai) in enumerate(zip(x_elec, q, eta_elec)):
for j, (xj, qj) in enumerate(zip(x_elyt, q_elyt)):
dx = xj - xi
fij_box = force_component(dx, np.abs(dx), qi * qj, etai)
fij_img = 4 * force_component(dx, image_distances[i][j], qi * qj, etai)
f_elec[i] -= fij_box + fij_img
f_elyt[j] += fij_box + fij_img
# electrolyte-electrolyte
for i, (xi, qi) in enumerate(zip(x_elyt, q_elyt)):
for j, (xj, qj) in list(enumerate(zip(x_elyt, q_elyt)))[i + 1 :]:
dx = xj - xi
fij_box = force_component(dx, np.abs(dx), qi * qj)
fij_img = 4 * force_component(dx, image_elyt_distances[i][j], qi * qj)
f_elyt[i] -= fij_img + fij_box
f_elyt[j] += fij_img + fij_box
# force units
assert np.abs(np.sum(f_elec) + np.sum(f_elyt)) < 1e-8
f_elec *= COULOMB
f_elyt *= COULOMB
# Virial
volume = a**2 * LZ
virial = 0.0
for x, f in [(x_elec, f_elec), (x_elyt, f_elyt)]:
virial += np.dot(x, f)
pressure = NKTV2P * virial / volume
with open(f"plate_cap{name}.csv", "w") as f:
f.write(
"length, energy / kcal/mol, q1 / e, q2 / e, inv11 / A, inv12 / A, b1 / e/A, b2 / e/A\n"
"length, energy / kcal/mol, q1 / e, q2 / e, inv11 / A, inv12 / A"
+ ", b1 / e/A, b2 / e/A, felec1 / kcal/mol/A, felec2 / kcal/mol/A"
+ ", felyt1 / kcal/mol/A, felyt2 / kcal/mol/A, press\n"
)
f.write(
", ".join(
@ -93,7 +164,14 @@ for name, eta_elec in [("", [2.0, 2.0]), ("_eta_mix", [0.5, 3.0])]:
f"{inv[0, 1]:.10f}",
f"{b[0]:.8f}",
f"{b[1]:.8f}",
f"{f_elec[0]:.5f}",
f"{f_elec[1]:.5f}",
f"{f_elyt[0]:.5f}",
f"{f_elyt[1]:.5f}",
f"{pressure:.2f}",
]
)
+ "\n"
)
time_end = time.perf_counter()
print(f"{time_end - time_start:0.4f} seconds")

6
examples/PACKAGES/electrode/madelung/settings.mod
View File

@ -19,4 +19,8 @@ compute qtop top reduce sum v_q
compute compute_pe all pe
variable vpe equal c_compute_pe
variable charge equal c_qtop
fix fxprint all print 1 "${vpe}, ${charge}" file "out.csv"
compute press all pressure NULL virial
variable p3 equal c_press[3]
fix fxprint all print 1 "${vpe}, ${charge}, ${p3}" file "out.csv"
dump dump_forces all custom 1 forces.lammpstrj id fx fy fz

4064
examples/PACKAGES/electrode/piston/data.piston.final
View File

File diff suppressed because it is too large Load Diff

6
examples/PACKAGES/electrode/piston/in.piston
View File

@ -42,7 +42,11 @@ fix fxforce_au gold setforce 0.0 0.0 0.0
# equilibrate z-coordinate of upper electrode while keeping the electrode rigid
fix fxforce_wa wall setforce 0.0 0.0 NULL
fix fxpressure wall aveforce 0 0 -0.005246 # atomspheric pressure: area/force->nktv2p
variable atm equal 1/68568.415 # 1/force->nktv2p
variable area equal (xhi-xlo)*(yhi-ylo)
variable wall_force equal -v_atm*v_area/count(wall)
print "Wall force per atom: ${wall_force}"
fix fxpressure wall aveforce 0 0 ${wall_force} # atomspheric pressure: area/force->nktv2p
fix fxdrag wall viscous 100
fix fxrigid wall rigid/nve single

97
examples/PACKAGES/electrode/piston/log.1Dec2022.piston.g++.1 → examples/PACKAGES/electrode/piston/log.22May2024.piston.g++.1
View File

@ -1,4 +1,5 @@
LAMMPS (3 Nov 2022)
LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-217-g1909233c69-modified)
using 1 OpenMP thread(s) per MPI task
# The intention is to find the average position of one wall at atmospheric
# pressure. The output is the wall position over time which can be used to
# find the average position for a run with fixed wall position.
@ -40,8 +41,8 @@ Finding 1-2 1-3 1-4 neighbors ...
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.011 seconds
special bonds CPU = 0.000 seconds
read_data CPU = 0.012 seconds
# ----------------- Settings Section -----------------
@ -77,7 +78,13 @@ fix fxforce_au gold setforce 0.0 0.0 0.0
# equilibrate z-coordinate of upper electrode while keeping the electrode rigid
fix fxforce_wa wall setforce 0.0 0.0 NULL
fix fxpressure wall aveforce 0 0 -0.005246 # atomspheric pressure: area/force->nktv2p
variable atm equal 1/68568.415 # 1/force->nktv2p
variable area equal (xhi-xlo)*(yhi-ylo)
variable wall_force equal -v_atm*v_area/count(wall)
print "Wall force per atom: ${wall_force}"
Wall force per atom: -0.000109285996244287
fix fxpressure wall aveforce 0 0 ${wall_force} # atomspheric pressure: area/force->nktv2p
fix fxpressure wall aveforce 0 0 -0.000109285996244287
fix fxdrag wall viscous 100
fix fxrigid wall rigid/nve single
1 rigid bodies with 48 atoms
@ -134,7 +141,7 @@ PPPM/electrode initialization ...
stencil order = 5
estimated absolute RMS force accuracy = 0.02930901
estimated relative force accuracy = 8.8263214e-05
using double precision MKL FFT
using double precision FFTW3
3d grid and FFT values/proc = 15884 6480
Generated 6 of 6 mixed pair_coeff terms from arithmetic mixing rule
Neighbor list info ...
@ -157,54 +164,54 @@ Neighbor list info ...
Per MPI rank memory allocation (min/avg/max) = 11.7 | 11.7 | 11.7 Mbytes
Step c_temp_mobile c_qwa c_qau v_top_wall
0 303.38967 -0.042963484 0.042963484 21.4018
5000 285.08828 -0.26105255 0.26105255 25.155629
10000 323.19176 -0.26264003 0.26264003 24.541676
15000 310.479 -0.27318148 0.27318148 23.141522
20000 295.18544 -0.11313444 0.11313444 23.828735
25000 295.38607 -0.25433086 0.25433086 23.673314
30000 288.0613 -0.30099901 0.30099901 23.438086
35000 278.5591 -0.15823576 0.15823576 24.311915
40000 303.95751 -0.19941381 0.19941381 23.69594
45000 279.026 -0.1659962 0.1659962 23.588604
50000 298.79278 -0.28866703 0.28866703 23.372508
55000 301.03353 -0.078370381 0.078370381 23.192985
60000 306.77965 -0.12807205 0.12807205 23.968574
65000 309.86008 -0.27162663 0.27162663 23.616704
70000 287.31116 -0.029751882 0.029751882 23.667495
75000 312.48654 -0.10759866 0.10759866 23.504105
80000 309.94267 -0.2558548 0.2558548 23.810576
85000 328.04389 -0.1575471 0.1575471 24.013437
90000 302.9806 -0.032002164 0.032002164 24.264432
95000 294.20804 -0.27797238 0.27797238 23.291758
100000 307.63019 -0.19047448 0.19047448 23.632147
Loop time of 530.844 on 1 procs for 100000 steps with 726 atoms
5000 311.85363 0.03543775 -0.03543775 24.79665
10000 285.91321 -0.16873703 0.16873703 23.103088
15000 295.39476 -0.44424612 0.44424612 23.767107
20000 296.12969 -0.14120993 0.14120993 23.96361
25000 306.59629 -0.29333182 0.29333182 23.884488
30000 297.98559 -0.10749684 0.10749684 23.73316
35000 297.98503 -0.11809975 0.11809975 23.984669
40000 300.26292 -0.32784184 0.32784184 23.462748
45000 295.68441 -0.25940165 0.25940165 23.516403
50000 315.12883 -0.36037614 0.36037614 23.627879
55000 290.55151 -0.0032838106 0.0032838106 23.684931
60000 316.4625 -0.17245368 0.17245368 24.126883
65000 296.79343 -0.054061851 0.054061851 23.695094
70000 305.99923 -0.11363801 0.11363801 23.55476
75000 297.40131 -0.27054153 0.27054153 23.928994
80000 306.54811 -0.25409719 0.25409719 23.869448
85000 303.95231 -0.17895561 0.17895561 23.658833
90000 313.43739 -0.059036514 0.059036514 23.36056
95000 290.3077 -0.31394478 0.31394478 23.885538
100000 297.5156 -0.30730083 0.30730083 23.511674
Loop time of 1586.06 on 1 procs for 100000 steps with 726 atoms
Performance: 32.552 ns/day, 0.737 hours/ns, 188.379 timesteps/s, 136.763 katom-step/s
100.0% CPU use with 1 MPI tasks x no OpenMP threads
Performance: 10.895 ns/day, 2.203 hours/ns, 63.049 timesteps/s, 45.774 katom-step/s
99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 190.47 | 190.47 | 190.47 | 0.0 | 35.88
Bond | 0.10754 | 0.10754 | 0.10754 | 0.0 | 0.02
Kspace | 73.179 | 73.179 | 73.179 | 0.0 | 13.79
Neigh | 24.209 | 24.209 | 24.209 | 0.0 | 4.56
Comm | 1.6857 | 1.6857 | 1.6857 | 0.0 | 0.32
Output | 0.0016861 | 0.0016861 | 0.0016861 | 0.0 | 0.00
Modify | 240.23 | 240.23 | 240.23 | 0.0 | 45.26
Other | | 0.9595 | | | 0.18
Pair | 460.91 | 460.91 | 460.91 | 0.0 | 29.06
Bond | 0.047873 | 0.047873 | 0.047873 | 0.0 | 0.00
Kspace | 341.4 | 341.4 | 341.4 | 0.0 | 21.53
Neigh | 52.868 | 52.868 | 52.868 | 0.0 | 3.33
Comm | 5.2321 | 5.2321 | 5.2321 | 0.0 | 0.33
Output | 0.00099102 | 0.00099102 | 0.00099102 | 0.0 | 0.00
Modify | 724.63 | 724.63 | 724.63 | 0.0 | 45.69
Other | | 0.9741 | | | 0.06
Nlocal: 726 ave 726 max 726 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 2335 ave 2335 max 2335 min
Nghost: 2336 ave 2336 max 2336 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 120271 ave 120271 max 120271 min
Neighs: 120321 ave 120321 max 120321 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 120271
Ave neighs/atom = 165.66253
Total # of neighbors = 120321
Ave neighs/atom = 165.7314
Ave special neighs/atom = 1.7355372
Neighbor list builds = 7722
Neighbor list builds = 7670
Dangerous builds = 0
write_data "data.piston.final"
System init for write_data ...
@ -213,11 +220,11 @@ PPPM/electrode initialization ...
G vector (1/distance) = 0.32814871
grid = 12 15 36
stencil order = 5
estimated absolute RMS force accuracy = 0.029311365
estimated relative force accuracy = 8.8270304e-05
using double precision MKL FFT
estimated absolute RMS force accuracy = 0.029311329
estimated relative force accuracy = 8.8270197e-05
using double precision FFTW3
3d grid and FFT values/proc = 15884 6480
Generated 6 of 6 mixed pair_coeff terms from arithmetic mixing rule
Average conjugate gradient steps: 1.981
Total wall time: 0:08:50
Total wall time: 0:26:26

109
examples/PACKAGES/electrode/piston/log.1Dec2022.piston.g++.4 → examples/PACKAGES/electrode/piston/log.22May2024.piston.g++.4
View File

@ -1,4 +1,5 @@
LAMMPS (3 Nov 2022)
LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-217-g1909233c69-modified)
using 1 OpenMP thread(s) per MPI task
# The intention is to find the average position of one wall at atmospheric
# pressure. The output is the wall position over time which can be used to
# find the average position for a run with fixed wall position.
@ -41,8 +42,8 @@ Finding 1-2 1-3 1-4 neighbors ...
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.017 seconds
special bonds CPU = 0.000 seconds
read_data CPU = 0.012 seconds
# ----------------- Settings Section -----------------
@ -66,7 +67,7 @@ Finding SHAKE clusters ...
0 = # of size 3 clusters
0 = # of size 4 clusters
210 = # of frozen angles
find clusters CPU = 0.002 seconds
find clusters CPU = 0.000 seconds
pair_modify mix arithmetic
# ----------------- Run Section -----------------
@ -78,7 +79,13 @@ fix fxforce_au gold setforce 0.0 0.0 0.0
# equilibrate z-coordinate of upper electrode while keeping the electrode rigid
fix fxforce_wa wall setforce 0.0 0.0 NULL
fix fxpressure wall aveforce 0 0 -0.005246 # atomspheric pressure: area/force->nktv2p
variable atm equal 1/68568.415 # 1/force->nktv2p
variable area equal (xhi-xlo)*(yhi-ylo)
variable wall_force equal -v_atm*v_area/count(wall)
print "Wall force per atom: ${wall_force}"
Wall force per atom: -0.000109285996244287
fix fxpressure wall aveforce 0 0 ${wall_force} # atomspheric pressure: area/force->nktv2p
fix fxpressure wall aveforce 0 0 -0.000109285996244287
fix fxdrag wall viscous 100
fix fxrigid wall rigid/nve single
1 rigid bodies with 48 atoms
@ -135,7 +142,7 @@ PPPM/electrode initialization ...
stencil order = 5
estimated absolute RMS force accuracy = 0.02930901
estimated relative force accuracy = 8.8263214e-05
using double precision MKL FFT
using double precision FFTW3
3d grid and FFT values/proc = 8512 2880
Generated 6 of 6 mixed pair_coeff terms from arithmetic mixing rule
Neighbor list info ...
@ -158,54 +165,54 @@ Neighbor list info ...
Per MPI rank memory allocation (min/avg/max) = 10.06 | 10.22 | 10.41 Mbytes
Step c_temp_mobile c_qwa c_qau v_top_wall
0 303.38967 -0.042963484 0.042963484 21.4018
5000 292.03027 -0.19040435 0.19040435 24.581338
10000 309.52764 -0.48308301 0.48308301 23.776985
15000 295.00243 -0.16591109 0.16591109 23.672038
20000 293.5536 -0.086669084 0.086669084 23.426455
25000 303.0079 -0.16488112 0.16488112 23.862966
30000 306.31463 -0.23192653 0.23192653 23.819882
35000 303.66268 -0.2317907 0.2317907 23.495344
40000 301.39435 -0.34661329 0.34661329 23.657835
45000 291.61205 -0.30539427 0.30539427 23.437303
50000 298.65319 -0.096107034 0.096107034 23.57809
55000 282.65069 -0.14943539 0.14943539 23.823728
60000 310.64182 -0.17418813 0.17418813 23.286959
65000 308.47141 -0.02075662 0.02075662 23.91313
70000 292.5186 -0.080163162 0.080163162 23.96283
75000 270.13928 -0.029528648 0.029528648 23.488972
80000 322.10914 0.030761045 -0.030761045 23.47592
85000 310.60347 -0.24069996 0.24069996 23.987091
90000 294.35695 -0.070458235 0.070458235 23.397929
95000 308.69043 -0.2652581 0.2652581 23.473813
100000 318.71883 0.024035956 -0.024035956 23.449863
Loop time of 590.232 on 4 procs for 100000 steps with 726 atoms
5000 291.6303 -0.1820085 0.1820085 24.641399
10000 299.42886 -0.19823095 0.19823095 23.820522
15000 288.23071 -0.065261869 0.065261869 23.360845
20000 299.4644 -0.042993777 0.042993777 23.987554
25000 304.26497 -0.15665293 0.15665293 23.729006
30000 292.29674 -0.25142779 0.25142779 23.960725
35000 295.57492 -0.01269228 0.01269228 23.445383
40000 303.38438 -0.13941727 0.13941727 23.517483
45000 302.211 -0.19589892 0.19589892 23.704043
50000 281.64939 -0.18057298 0.18057298 23.542137
55000 274.90565 -0.15453379 0.15453379 23.734347
60000 290.70459 -0.27977436 0.27977436 23.835365
65000 293.42241 -0.2454241 0.2454241 23.59269
70000 295.20229 -0.041314995 0.041314995 23.73856
75000 297.79519 -0.11231755 0.11231755 23.57262
80000 285.17858 -0.070796508 0.070796508 23.817135
85000 311.71609 -0.068920177 0.068920177 23.861127
90000 287.80446 -0.19183387 0.19183387 23.369393
95000 309.43345 -0.15238671 0.15238671 23.597792
100000 294.12422 -0.14284353 0.14284353 23.526286
Loop time of 876.546 on 4 procs for 100000 steps with 726 atoms
Performance: 29.277 ns/day, 0.820 hours/ns, 169.425 timesteps/s, 123.003 katom-step/s
72.5% CPU use with 4 MPI tasks x no OpenMP threads
Performance: 19.714 ns/day, 1.217 hours/ns, 114.084 timesteps/s, 82.825 katom-step/s
98.6% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 57.391 | 75.867 | 96.292 | 212.1 | 12.85
Bond | 0.10177 | 0.11042 | 0.12415 | 2.7 | 0.02
Kspace | 102.79 | 123.16 | 141.5 | 165.7 | 20.87
Neigh | 12.808 | 12.895 | 12.982 | 2.3 | 2.18
Comm | 18.885 | 19.973 | 21.064 | 24.0 | 3.38
Output | 0.0022573 | 0.0022749 | 0.0023225 | 0.1 | 0.00
Modify | 355.89 | 356.74 | 357.61 | 4.2 | 60.44
Other | | 1.478 | | | 0.25
Pair | 123.63 | 171.23 | 215.73 | 336.6 | 19.53
Bond | 0.068261 | 0.075883 | 0.081822 | 1.9 | 0.01
Kspace | 187.59 | 231.71 | 279.01 | 287.1 | 26.43
Neigh | 29.28 | 29.462 | 29.637 | 2.5 | 3.36
Comm | 12.544 | 13.731 | 14.929 | 29.1 | 1.57
Output | 0.0010182 | 0.0014585 | 0.0016071 | 0.7 | 0.00
Modify | 428.74 | 429.25 | 429.74 | 2.3 | 48.97
Other | | 1.092 | | | 0.12
Nlocal: 181.5 ave 207 max 169 min
Histogram: 2 0 1 0 0 0 0 0 0 1
Nghost: 1961.5 ave 1984 max 1926 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Neighs: 30051 ave 41646 max 20775 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Nlocal: 181.5 ave 195 max 166 min
Histogram: 1 1 0 0 0 0 0 0 0 2
Nghost: 1955.5 ave 1978 max 1931 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Neighs: 30343 ave 39847 max 20428 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 120204
Ave neighs/atom = 165.57025
Total # of neighbors = 121372
Ave neighs/atom = 167.17906
Ave special neighs/atom = 1.7355372
Neighbor list builds = 7663
Neighbor list builds = 7698
Dangerous builds = 0
write_data "data.piston.final"
System init for write_data ...
@ -214,11 +221,11 @@ PPPM/electrode initialization ...
G vector (1/distance) = 0.32814871
grid = 12 15 36
stencil order = 5
estimated absolute RMS force accuracy = 0.029311028
estimated relative force accuracy = 8.8269289e-05
using double precision MKL FFT
estimated absolute RMS force accuracy = 0.029310954
estimated relative force accuracy = 8.8269069e-05
using double precision FFTW3
3d grid and FFT values/proc = 8512 2880
Generated 6 of 6 mixed pair_coeff terms from arithmetic mixing rule
Average conjugate gradient steps: 1.982
Total wall time: 0:09:50
Average conjugate gradient steps: 1.981
Total wall time: 0:14:36

75
examples/rheo/balloon/in.rheo.balloon Normal file
View File

@ -0,0 +1,75 @@
# ------ 2D water balloon ------ #
dimension 2
units lj
atom_style hybrid rheo bond
boundary m m p
comm_modify vel yes
newton off
region box block -40 40 0 80 -0.01 0.01 units box
create_box 1 box bond/types 1 extra/bond/per/atom 15 extra/special/per/atom 50
region fluid sphere -10 40 0 30 units box side in
lattice hex 1.0
create_atoms 1 region fluid
region shell sphere -10 40 0 27 units box side out
group shell region shell
set group shell rheo/status 1
set group all vx 0.005 vy -0.04
# ------ Model parameters ------#
variable cut equal 3.0
variable n equal 1.0
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable zeta equal 0.05
variable kappa equal 0.01*${rho0}/${mp}
variable dt_max equal 0.1*${cut}/${cs}/3
variable eta equal 0.05
variable Cv equal 1.0
variable L equal 1.0
variable Tf equal 1.0
mass * ${mp}
timestep 0.1
pair_style hybrid/overlay rheo ${cut} artificial/visc ${zeta} rheo/solid
pair_coeff * * rheo
pair_coeff * * rheo/solid 1.0 1.0 1.0
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
create_bonds many shell shell 1 0 1.5
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
bond_style bpm/spring
bond_coeff 1 1.0 1.0 1.0
# A lower critical strain allows the balloon to pop
#bond_coeff 1 1.0 0.05 1.0
# ------ Drop balloon ------#
fix 1 all rheo ${cut} quintic 0 &
shift &
surface/detection coordination 22 8
fix 2 all rheo/viscosity * constant ${eta}
fix 3 all rheo/pressure * linear
fix 4 all wall/harmonic ylo EDGE 2.0 1.0 1.0
fix 5 all enforce2d
compute rho all rheo/property/atom rho
compute phase all rheo/property/atom phase
compute nbond all nbond/atom
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press atoms
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_phase c_nbond c_rho
run 30000

382
examples/rheo/balloon/log.17Apr2024.balloon.g++.4 Normal file
View File

@ -0,0 +1,382 @@
LAMMPS (17 Apr 2024 - Development - patch_5May2020-18508-g3c0eaf6870-modified)
# ------ 2D water balloon ------ #
dimension 2
units lj
atom_style hybrid rheo bond
boundary m m p
comm_modify vel yes
newton off
region box block -40 40 0 80 -0.01 0.01 units box
create_box 1 box bond/types 1 extra/bond/per/atom 15 extra/special/per/atom 50
Created orthogonal box = (-40 0 -0.01) to (40 80 0.01)
2 by 2 by 1 MPI processor grid
region fluid sphere -10 40 0 30 units box side in
lattice hex 1.0
Lattice spacing in x,y,z = 1.0745699 1.8612097 1.0745699
create_atoms 1 region fluid
Created 2830 atoms
using lattice units in orthogonal box = (-40 0 -0.01) to (40 80 0.01)
create_atoms CPU = 0.001 seconds
region shell sphere -10 40 0 27 units box side out
group shell region shell
544 atoms in group shell
set group shell rheo/status 1
Setting atom values ...
544 settings made for rheo/status
set group all vx 0.005 vy -0.04
Setting atom values ...
2830 settings made for vx
2830 settings made for vy
# ------ Model parameters ------#
variable cut equal 3.0
variable n equal 1.0
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable mp equal 1/${n}
variable mp equal 1/1
variable zeta equal 0.05
variable kappa equal 0.01*${rho0}/${mp}
variable kappa equal 0.01*1/${mp}
variable kappa equal 0.01*1/1
variable dt_max equal 0.1*${cut}/${cs}/3
variable dt_max equal 0.1*3/${cs}/3
variable dt_max equal 0.1*3/1/3
variable eta equal 0.05
variable Cv equal 1.0
variable L equal 1.0
variable Tf equal 1.0
mass * ${mp}
mass * 1
timestep 0.1
pair_style hybrid/overlay rheo ${cut} artificial/visc ${zeta} rheo/solid
pair_style hybrid/overlay rheo 3 artificial/visc ${zeta} rheo/solid
pair_style hybrid/overlay rheo 3 artificial/visc 0.05 rheo/solid
pair_coeff * * rheo
pair_coeff * * rheo/solid 1.0 1.0 1.0
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
0 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.000 seconds
create_bonds many shell shell 1 0 1.5
Generated 0 of 0 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 = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 49 49 1
3 neighbor lists, perpetual/occasional/extra = 2 1 0
(1) command create_bonds, occasional
attributes: full, newton off
pair build: full/bin
stencil: full/bin/2d
bin: standard
(2) pair rheo, perpetual
attributes: half, newton off
pair build: half/bin/newtoff
stencil: full/bin/2d
bin: standard
(3) pair rheo/solid, perpetual, trim from (2)
attributes: half, newton off, cut 1.3
pair build: trim
stencil: none
bin: none
Added 1263 bonds, new total = 1263
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
6 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.000 seconds
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
bond_style bpm/spring
bond_coeff 1 1.0 1.0 1.0
# A lower critical strain allows the balloon to pop
#bond_coeff 1 1.0 0.05 1.0
# ------ Drop balloon ------#
fix 1 all rheo ${cut} quintic 0 shift surface/detection coordination 22 8
fix 1 all rheo 3 quintic 0 shift surface/detection coordination 22 8
fix 2 all rheo/viscosity * constant ${eta}
fix 2 all rheo/viscosity * constant 0.05
fix 3 all rheo/pressure * linear
fix 4 all wall/harmonic ylo EDGE 2.0 1.0 1.0
fix 5 all enforce2d
compute rho all rheo/property/atom rho
compute phase all rheo/property/atom phase
compute nbond all nbond/atom
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press atoms
dump 1 all custom 200 atomDump id type x y vx vy fx fy c_phase c_nbond c_rho
run 30000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- BPM bond style: doi:10.1039/D3SM01373A
@Article{Clemmer2024,
author = {Clemmer, Joel T. and Monti, Joseph M. and Lechman, Jeremy B.},
title = {A soft departure from jamming: the compaction of deformable
granular matter under high pressures},
journal = {Soft Matter},
year = 2024,
volume = 20,
number = 8,
pages = {1702--1718}
}
- @article{PalermoInPrep,
journal = {in prep},
title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},
year = {2024},
author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 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 = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 49 49 1
6 neighbor lists, perpetual/occasional/extra = 6 0 0
(1) pair rheo, perpetual, half/full from (3)
attributes: half, newton off
pair build: halffull/newtoff
stencil: none
bin: none
(2) pair rheo/solid, perpetual, trim from (1)
attributes: half, newton off, cut 1.3
pair build: trim
stencil: none
bin: none
(3) compute RHEO/KERNEL, perpetual
attributes: full, newton off
pair build: full/bin
stencil: full/bin/2d
bin: standard
(4) compute RHEO/GRAD, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(5) compute RHEO/VSHIFT, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(6) compute RHEO/SURFACE, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 17.63 | 17.64 | 17.65 Mbytes
Step Time KinEng Press Atoms
0 0 0.0008125 0.00035927734 2830
200 20 0.0008125 0.00035927734 2830
400 40 0.0008125 0.00035927734 2830
600 60 0.0008125 0.00035927734 2830
800 80 0.0008125 0.00035927734 2830
1000 100 0.0008125 0.00035927734 2830
1200 120 0.0008125 0.00035927734 2830
1400 140 0.0008125 0.00035927734 2830
1600 160 0.0008125 0.00035927734 2830
1800 180 0.0008125 0.00035927734 2830
2000 200 0.0008125 0.00035927734 2830
2200 220 0.0008125 0.00035927734 2830
2400 240 0.00079033569 0.00043037861 2830
2600 260 0.0007549229 0.00045188383 2830
2800 280 0.00072808836 0.00031695003 2830
3000 300 0.0007017958 1.6121754e-05 2830
3200 320 0.00067479047 -0.00015725514 2830
3400 340 0.00064762254 -0.00023361314 2830
3600 360 0.00061960255 -0.00033837679 2830
3800 380 0.0005857206 -0.00051770716 2830
4000 400 0.00055061733 -0.00070309251 2830
4200 420 0.00051884719 -0.0008247795 2830
4400 440 0.00049022236 -0.00099918413 2830
4600 460 0.00046060011 -0.0010923159 2830
4800 480 0.00042900173 -0.0011524571 2830
5000 500 0.00039751503 -0.0012586358 2830
5200 520 0.00036620054 -0.0013973543 2830
5400 540 0.00033130023 -0.0015185231 2830
5600 560 0.00030565892 -0.0016159836 2830
5800 580 0.00028209836 -0.0016925198 2830
6000 600 0.00024695044 -0.0017796892 2830
6200 620 0.00021190635 -0.0018706272 2830
6400 640 0.0001947093 -0.0019146643 2830
6600 660 0.00018903936 -0.0019146199 2830
6800 680 0.00017753371 -0.0019390155 2830
7000 700 0.00015170593 -0.0020247472 2830
7200 720 0.00011509692 -0.0021222209 2830
7400 740 7.9861785e-05 -0.0022033181 2830
7600 760 6.1350463e-05 -0.0022511971 2830
7800 780 6.5269523e-05 -0.0022222806 2830
8000 800 8.5709569e-05 -0.0021089664 2830
8200 820 0.00011746348 -0.0019351493 2830
8400 840 0.00015698134 -0.0017079928 2830
8600 860 0.00019758065 -0.0014618965 2830
8800 880 0.00023338199 -0.0012365832 2830
9000 900 0.00026282353 -0.0010348527 2830
9200 920 0.00028604776 -0.00085287884 2830
9400 940 0.00030388767 -0.000681122 2830
9600 960 0.000317589 -0.00052203521 2830
9800 980 0.00032716728 -0.00037501187 2830
10000 1000 0.00033270692 -0.00025576132 2830
10200 1020 0.00033485986 -0.00016554207 2830
10400 1040 0.00033476763 -9.8525417e-05 2830
10600 1060 0.00033351922 -5.1166347e-05 2830
10800 1080 0.00033161645 -2.0773965e-05 2830
11000 1100 0.00032913022 2.2384421e-07 2830
11200 1120 0.00032618376 1.2304773e-05 2830
11400 1140 0.00032310409 1.3725982e-05 2830
11600 1160 0.0003202128 9.0431945e-06 2830
11800 1180 0.00031760386 -5.3537879e-07 2830
12000 1200 0.00031518884 -1.331708e-05 2830
12200 1220 0.00031283958 -3.0838612e-05 2830
12400 1240 0.0003104901 -5.0038548e-05 2830
12600 1260 0.00030811597 -6.9699925e-05 2830
12800 1280 0.00030555782 -8.9972287e-05 2830
13000 1300 0.00030256671 -0.00011712941 2830
13200 1320 0.00029907961 -0.00015495826 2830
13400 1340 0.00029504656 -0.00020292633 2830
13600 1360 0.0002905184 -0.00024892421 2830
13800 1380 0.00028564542 -0.000295085 2830
14000 1400 0.00028073246 -0.00034571956 2830
14200 1420 0.00027611457 -0.00039341977 2830
14400 1440 0.00027217382 -0.0004281012 2830
14600 1460 0.00026919129 -0.00045342545 2830
14800 1480 0.00026727674 -0.00047323419 2830
15000 1500 0.0002663482 -0.00048423944 2830
15200 1520 0.00026616663 -0.0004816085 2830
15400 1540 0.00026634862 -0.00047573486 2830
15600 1560 0.0002664314 -0.00046803192 2830
15800 1580 0.00026603348 -0.00045753668 2830
16000 1600 0.00026511015 -0.00044676105 2830
16200 1620 0.00026373403 -0.00044075794 2830
16400 1640 0.00026217342 -0.00043684036 2830
16600 1660 0.0002607038 -0.00042774771 2830
16800 1680 0.00025951097 -0.00041603026 2830
17000 1700 0.00025869088 -0.00040302996 2830
17200 1720 0.00025825588 -0.00038415247 2830
17400 1740 0.00025818373 -0.00035742127 2830
17600 1760 0.00025843381 -0.00032854722 2830
17800 1780 0.00025897836 -0.00029821183 2830
18000 1800 0.00025981472 -0.00026108907 2830
18200 1820 0.00026095775 -0.00021731058 2830
18400 1840 0.00026239688 -0.00017030825 2830
18600 1860 0.00026404432 -0.00011868778 2830
18800 1880 0.00026574247 -5.9556286e-05 2830
19000 1900 0.00026729563 2.3014881e-06 2830
19200 1920 0.00026852418 6.2100169e-05 2830
19400 1940 0.00026929086 0.00012090325 2830
19600 1960 0.0002695407 0.00017904223 2830
19800 1980 0.00026929677 0.00023112254 2830
20000 2000 0.00026863577 0.0002756697 2830
20200 2020 0.00026765699 0.0003158399 2830
20400 2040 0.00026646841 0.00035200747 2830
20600 2060 0.00026516938 0.00038018442 2830
20800 2080 0.00026383495 0.00040179111 2830
21000 2100 0.00026252489 0.00042030921 2830
21200 2120 0.00026128616 0.00043466976 2830
21400 2140 0.00026014896 0.00044221445 2830
21600 2160 0.00025912325 0.00044531883 2830
21800 2180 0.00025821515 0.00044661709 2830
22000 2200 0.00025742576 0.00044409089 2830
22200 2220 0.00025674938 0.00043634999 2830
22400 2240 0.00025617111 0.00042630344 2830
22600 2260 0.0002556791 0.00041561603 2830
22800 2280 0.00025525963 0.00040166735 2830
23000 2300 0.00025489538 0.00038430419 2830
23200 2320 0.00025456861 0.0003669402 2830
23400 2340 0.00025426747 0.00034972373 2830
23600 2360 0.00025398353 0.0003302242 2830
23800 2380 0.00025370842 0.00030993088 2830
24000 2400 0.00025344084 0.00029143258 2830
24200 2420 0.00025318683 0.00027421708 2830
24400 2440 0.0002529591 0.00025603123 2830
24600 2460 0.0002527713 0.00023950245 2830
24800 2480 0.00025264228 0.00022644812 2830
25000 2500 0.00025259021 0.00021540748 2830
25200 2520 0.00025262892 0.00020544201 2830
25400 2540 0.00025276229 0.00019845807 2830
25600 2560 0.0002529876 0.00019449958 2830
25800 2580 0.00025329374 0.00019082606 2830
26000 2600 0.00025366066 0.00018700064 2830
26200 2620 0.00025406164 0.00018426061 2830
26400 2640 0.00025446737 0.00018098339 2830
26600 2660 0.00025484714 0.00017471869 2830
26800 2680 0.00025516604 0.00016565557 2830
27000 2700 0.00025538911 0.00015493626 2830
27200 2720 0.00025548177 0.00014075592 2830
27400 2740 0.00025541168 0.00012205573 2830
27600 2760 0.00025514889 0.00010039772 2830
27800 2780 0.00025467547 7.7069215e-05 2830
28000 2800 0.0002539915 5.1158042e-05 2830
28200 2820 0.00025312083 2.3468384e-05 2830
28400 2840 0.00025211323 -3.2184465e-06 2830
28600 2860 0.00025104366 -2.7726301e-05 2830
28800 2880 0.00025000263 -5.0202987e-05 2830
29000 2900 0.00024907814 -6.9244776e-05 2830
29200 2920 0.00024833815 -8.2874516e-05 2830
29400 2940 0.0002478155 -9.1854992e-05 2830
29600 2960 0.00024750313 -9.766055e-05 2830
29800 2980 0.00024735538 -9.9681291e-05 2830
30000 3000 0.00024730191 -9.818759e-05 2830
Loop time of 177.982 on 4 procs for 30000 steps with 2830 atoms
Performance: 1456330.235 tau/day, 168.557 timesteps/s, 477.016 katom-step/s
99.7% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 22.913 | 27.061 | 34.594 | 87.2 | 15.20
Bond | 0.22386 | 0.26159 | 0.30792 | 6.0 | 0.15
Neigh | 0.84412 | 0.84509 | 0.8462 | 0.1 | 0.47
Comm | 0.50015 | 0.55579 | 0.60346 | 5.2 | 0.31
Output | 0.65854 | 0.69412 | 0.72473 | 2.8 | 0.39
Modify | 133.13 | 136 | 137.38 | 14.5 | 76.41
Other | | 12.57 | | | 7.06
Nlocal: 707.5 ave 1576 max 53 min
Histogram: 2 0 0 0 0 0 1 0 0 1
Nghost: 164.75 ave 239 max 94 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Neighs: 12307.8 ave 27380 max 983 min
Histogram: 2 0 0 0 0 0 1 0 0 1
FullNghs: 23517 ave 53040 max 1502 min
Histogram: 2 0 0 0 0 0 1 0 0 1
Total # of neighbors = 94068
Ave neighs/atom = 33.239576
Ave special neighs/atom = 0.89257951
Neighbor list builds = 783
Dangerous builds = 0
Total wall time: 0:02:58

76
examples/rheo/dam-break/in.rheo.dam.break Normal file
View File

@ -0,0 +1,76 @@
# ------ 2D dam break ------ #
dimension 2
units lj
atom_style rheo
boundary f s p
comm_modify vel yes
newton off
# ------ Create simulation box ------ #
variable n equal 1.0
variable cut equal 2.2
variable dx equal 3.0
region box block -1 150 -1 80 -0.1 0.1 units box
create_box 2 box
lattice hex ${n}
region fluid block $(xlo+v_dx+1.0) $(xlo+40.0) $(ylo+v_dx+1.0) $(yhi-20.0) EDGE EDGE units box
region walls1 block $(xlo+v_dx) $(xhi-v_dx) $(ylo+v_dx) $(yhi-v_dx) EDGE EDGE side out units box
region walls2 block EDGE EDGE EDGE $(yhi-v_dx) EDGE EDGE side in units box
region walls intersect 2 walls1 walls2
create_atoms 1 region fluid
create_atoms 2 region walls
group fluid type 1
group rig type 2
# ------ Model parameters ------ #
variable rho0 equal 1.0
variable mp equal ${rho0}/${n}
variable cs equal 1.0
variable zeta equal 0.1
variable dt_max equal 0.1*${cut}/${cs}/3
variable eta equal 0.1
variable Dr equal 0.1
mass 1 ${mp}
mass 2 $(2*v_mp)
set group all rheo/rho ${rho0}
set group all rheo/status 0
set group rig rheo/status 1
timestep ${dt_max}
pair_style rheo ${cut} artificial/visc ${zeta} #rho/damp ${Dr}
pair_coeff * *
# ------ Fixes & computes ------ #
fix 1 all rheo ${cut} quintic 10 &
surface/detection coordination 22 8 &
rho/sum
fix 2 all rheo/viscosity * constant ${eta}
fix 3 all rheo/pressure * linear
fix 4 all gravity 1e-3 vector 0 -1 0
fix 5 rig setforce 0.0 0.0 0.0
fix 6 all enforce2d
compute rho all rheo/property/atom rho
compute p all rheo/property/atom pressure
compute surf all rheo/property/atom surface
compute sn all rheo/property/atom surface/n/x surface/n/y
# ------ Output & Run ------ #
thermo 20
thermo_style custom step time ke press
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_rho c_surf c_p c_sn[*]
run 30000

1694
examples/rheo/dam-break/log.17Apr2024.dam.g++.4 Normal file
View File

File diff suppressed because it is too large Load Diff

82
examples/rheo/ice-cubes/in.rheo.ice.cubes Normal file
View File

@ -0,0 +1,82 @@
# ------ 2D Ice Cube Pour ------ #
dimension 2
units lj
atom_style hybrid rheo/thermal bond
boundary m m p
comm_modify vel yes
newton off
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
region box block -25 25 0 100 -0.01 0.01 units box
create_box 1 box bond/types 1 extra/bond/per/atom 15 extra/special/per/atom 50
region fluid block $(xlo+1) $(xhi-1) $(ylo+1) $(ylo+30) EDGE EDGE units box
lattice sq 1.0
create_atoms 1 region fluid
set group all sph/e 8.0
# ------ Model parameters ------#
variable cut equal 3.0
variable n equal 1.0
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable zeta equal 0.05
variable kappa equal 0.01*${rho0}/${mp}
variable dt_max equal 0.1*${cut}/${cs}/3
variable eta equal 0.05
variable Cv equal 1.0
variable L equal 1.0
variable Tf equal 1.0
mass * ${mp}
timestep 0.1
pair_style hybrid/overlay rheo ${cut} artificial/visc ${zeta} rheo/solid
pair_coeff * * rheo
pair_coeff * * rheo/solid 1.0 1.0 1.0
bond_style bpm/spring
bond_coeff 1 1.0 1.0 1.0
# ------ Pour particles ------#
molecule my_mol "square.mol"
# Wall region extends far enough in z to avoid contact
region wall block EDGE EDGE EDGE EDGE -5 5 side in open 4 units box
region drop block -16 16 70 90 EDGE EDGE side in units box
fix 1 all rheo ${cut} quintic 0 &
thermal &
shift &
surface/detection coordination 22 8
fix 2 all rheo/viscosity * constant ${eta}
fix 3 all rheo/pressure * linear
fix 4 all rheo/thermal conductivity * constant ${kappa} &
specific/heat * constant ${Cv} &
Tfreeze * constant ${Tf} &
latent/heat * constant ${L} &
react 1.5 1
fix 5 all wall/region wall harmonic 1.0 1.0 1.0
fix 6 all gravity 5e-4 vector 0 -1 0
fix 7 all deposit 8 0 1000 37241459 mol my_mol region drop near 2.0 vy -0.02 -0.02
fix 8 all enforce2d
compute rho all rheo/property/atom rho
compute phase all rheo/property/atom phase
compute temp all rheo/property/atom temperature
compute eng all rheo/property/atom energy
compute nbond all nbond/atom
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press atoms
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_phase c_temp c_eng c_nbond c_rho
run 30000

379
examples/rheo/ice-cubes/log.17Apr2024.ice.g++.4 Normal file
View File

@ -0,0 +1,379 @@
LAMMPS (17 Apr 2024 - Development - patch_5May2020-18508-g3c0eaf6870-modified)
# ------ 2D Ice Cube Pour ------ #
dimension 2
units lj
atom_style hybrid rheo/thermal bond
boundary m m p
comm_modify vel yes
newton off
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
region box block -25 25 0 100 -0.01 0.01 units box
create_box 1 box bond/types 1 extra/bond/per/atom 15 extra/special/per/atom 50
Created orthogonal box = (-25 0 -0.01) to (25 100 0.01)
2 by 2 by 1 MPI processor grid
region fluid block $(xlo+1) $(xhi-1) $(ylo+1) $(ylo+30) EDGE EDGE units box
region fluid block -24 $(xhi-1) $(ylo+1) $(ylo+30) EDGE EDGE units box
region fluid block -24 24 $(ylo+1) $(ylo+30) EDGE EDGE units box
region fluid block -24 24 1 $(ylo+30) EDGE EDGE units box
region fluid block -24 24 1 30 EDGE EDGE units box
lattice sq 1.0
Lattice spacing in x,y,z = 1 1 1
create_atoms 1 region fluid
Created 1470 atoms
using lattice units in orthogonal box = (-25 0 -0.01) to (25 100 0.01)
create_atoms CPU = 0.001 seconds
set group all sph/e 8.0
Setting atom values ...
1470 settings made for sph/e
# ------ Model parameters ------#
variable cut equal 3.0
variable n equal 1.0
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable mp equal 1/${n}
variable mp equal 1/1
variable zeta equal 0.05
variable kappa equal 0.01*${rho0}/${mp}
variable kappa equal 0.01*1/${mp}
variable kappa equal 0.01*1/1
variable dt_max equal 0.1*${cut}/${cs}/3
variable dt_max equal 0.1*3/${cs}/3
variable dt_max equal 0.1*3/1/3
variable eta equal 0.05
variable Cv equal 1.0
variable L equal 1.0
variable Tf equal 1.0
mass * ${mp}
mass * 1
timestep 0.1
pair_style hybrid/overlay rheo ${cut} artificial/visc ${zeta} rheo/solid
pair_style hybrid/overlay rheo 3 artificial/visc ${zeta} rheo/solid
pair_style hybrid/overlay rheo 3 artificial/visc 0.05 rheo/solid
pair_coeff * * rheo
pair_coeff * * rheo/solid 1.0 1.0 1.0
bond_style bpm/spring
bond_coeff 1 1.0 1.0 1.0
# ------ Pour particles ------#
molecule my_mol "square.mol"
Read molecule template my_mol:
#Made with create_mol.py
1 molecules
0 fragments
100 atoms with max type 1
342 bonds with max type 1
0 angles with max type 0
0 dihedrals with max type 0
0 impropers with max type 0
# Wall region extends far enough in z to avoid contact
region wall block EDGE EDGE EDGE EDGE -5 5 side in open 4 units box
region drop block -16 16 70 90 EDGE EDGE side in units box
fix 1 all rheo ${cut} quintic 0 thermal shift surface/detection coordination 22 8
fix 1 all rheo 3 quintic 0 thermal shift surface/detection coordination 22 8
fix 2 all rheo/viscosity * constant ${eta}
fix 2 all rheo/viscosity * constant 0.05
fix 3 all rheo/pressure * linear
fix 4 all rheo/thermal conductivity * constant ${kappa} specific/heat * constant ${Cv} Tfreeze * constant ${Tf} latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.01 specific/heat * constant ${Cv} Tfreeze * constant ${Tf} latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.01 specific/heat * constant 1 Tfreeze * constant ${Tf} latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.01 specific/heat * constant 1 Tfreeze * constant 1 latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.01 specific/heat * constant 1 Tfreeze * constant 1 latent/heat * constant 1 react 1.5 1
fix 5 all wall/region wall harmonic 1.0 1.0 1.0
fix 6 all gravity 5e-4 vector 0 -1 0
fix 7 all deposit 8 0 1000 37241459 mol my_mol region drop near 2.0 vy -0.02 -0.02
WARNING: Molecule attributes do not match system attributes (../molecule.cpp:1881)
fix 8 all enforce2d
compute rho all rheo/property/atom rho
compute phase all rheo/property/atom phase
compute temp all rheo/property/atom temperature
compute eng all rheo/property/atom energy
compute nbond all nbond/atom
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press atoms
dump 1 all custom 200 atomDump id type x y vx vy fx fy c_phase c_temp c_eng c_nbond c_rho
run 30000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- BPM bond style: doi:10.1039/D3SM01373A
@Article{Clemmer2024,
author = {Clemmer, Joel T. and Monti, Joseph M. and Lechman, Jeremy B.},
title = {A soft departure from jamming: the compaction of deformable
granular matter under high pressures},
journal = {Soft Matter},
year = 2024,
volume = 20,
number = 8,
pages = {1702--1718}
}
- @article{PalermoInPrep,
journal = {in prep},
title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},
year = {2024},
author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},
}
- @article{ApplMathModel.130.310,
title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},
journal = {Applied Mathematical Modelling},
volume = {130},
pages = {310-326},
year = {2024},
issn = {0307-904X},
doi = {https://doi.org/10.1016/j.apm.2024.02.027},
author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 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 = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 31 61 1
7 neighbor lists, perpetual/occasional/extra = 6 1 0
(1) pair rheo, perpetual, half/full from (3)
attributes: half, newton off
pair build: halffull/newtoff
stencil: none
bin: none
(2) pair rheo/solid, perpetual, trim from (4)
attributes: half, newton off, cut 1.3
pair build: trim
stencil: none
bin: none
(3) compute RHEO/KERNEL, perpetual
attributes: full, newton off
pair build: full/bin
stencil: full/bin/2d
bin: standard
(4) compute RHEO/GRAD, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(5) compute RHEO/VSHIFT, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(6) compute RHEO/SURFACE, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(7) fix rheo/thermal, occasional, trim from (4)
attributes: half, newton off, cut 3
pair build: trim
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 15.53 | 15.61 | 15.69 Mbytes
Step Time KinEng Press Atoms
0 0 0 0 1470
200 20 5.6002982e-05 3.4434234e-05 1570
400 40 8.2173099e-05 8.6171768e-05 1570
600 60 8.019018e-05 0.00010750355 1570
800 80 0.00013866953 0.00010265608 1570
1000 100 0.00018965028 8.1985605e-05 1570
1200 120 0.00022033242 7.4736443e-05 1670
1400 140 0.00030767062 0.00011264333 1670
1600 160 0.00040770127 0.00018779992 1670
1800 180 0.00047476332 0.00023153009 1670
2000 200 0.00059116774 0.00027200445 1670
2200 220 0.0007151733 0.0002919963 1770
2400 240 0.00083392135 0.00029757889 1770
2600 260 0.00099653466 0.00036547269 1770
2800 280 0.0011964069 0.00045983458 1770
3000 300 0.0013716953 0.00055013647 1770
3200 320 0.0015174096 0.00064203572 1870
3400 340 0.0016539743 0.00086671622 1870
3600 360 0.0015887858 0.00066353749 1870
3800 380 0.0016451684 0.00070551483 1870
4000 400 0.0017330971 0.00080722283 1870
4200 420 0.001812193 0.00073573903 1970
4400 440 0.001755871 0.0010621909 1970
4600 460 0.0016190772 0.00072913706 1970
4800 480 0.0015741931 0.00073524088 1970
5000 500 0.0016488815 0.00088684275 1970
5200 520 0.0017213288 0.00077042378 2070
5400 540 0.0018509598 0.0010219434 2070
5600 560 0.0020251064 0.00083182483 2070
5800 580 0.0022473255 0.00095076144 2070
6000 600 0.0024843519 0.0011247014 2070
6200 620 0.0022282321 0.0018105932 2170
6400 640 0.0020289063 0.0014158497 2170
6600 660 0.002145241 0.0011359383 2170
6800 680 0.0024313937 0.0016475504 2170
7000 700 0.0021000599 0.0020983745 2170
7200 720 0.0019137235 0.0010439152 2270
7400 740 0.0018801367 0.00095436448 2270
7600 760 0.0017979449 0.0011184039 2270
7800 780 0.0018005205 0.0009243205 2270
8000 800 0.0017827073 0.0013671228 2270
8200 820 0.0018387108 0.0015426012 2270
8400 840 0.0016000788 0.0016751514 2270
8600 860 0.0013954964 0.0016884335 2270
8800 880 0.0013283728 0.0012399398 2270
9000 900 0.001389385 0.0012968496 2270
9200 920 0.0012295438 0.0012995821 2270
9400 940 0.0010522655 0.00082245528 2270
9600 960 0.00097085496 0.00053833131 2270
9800 980 0.0009398987 0.00063467387 2270
10000 1000 0.00092710392 0.00059494446 2270
10200 1020 0.00095545471 0.00074560644 2270
10400 1040 0.0009645841 0.00085429807 2270
10600 1060 0.00064037148 0.0017222246 2270
10800 1080 0.00046790978 0.00088204234 2270
11000 1100 0.00030106229 0.00074950209 2270
11200 1120 0.00027746016 0.00052831745 2270
11400 1140 0.0002533348 0.0006272715 2270
11600 1160 0.00021825085 0.00029691552 2270
11800 1180 0.0001451308 0.00015037478 2270
12000 1200 0.0001314823 0.00017227174 2270
12200 1220 0.00013693632 0.00017791384 2270
12400 1240 0.00014987347 0.0002286677 2270
12600 1260 0.00015092598 0.0003698436 2270
12800 1280 0.0001291653 0.00047229532 2270
13000 1300 0.00011949988 0.00049560375 2270
13200 1320 0.00011694665 0.00057542084 2270
13400 1340 9.6164519e-05 0.00062714755 2270
13600 1360 8.4517591e-05 0.00044156913 2270
13800 1380 0.00019140516 0.0003264745 2270
14000 1400 0.00013868599 0.00037753497 2270
14200 1420 9.3701636e-05 0.00031517848 2270
14400 1440 6.7389077e-05 0.0002946861 2270
14600 1460 5.3640086e-05 0.00026650711 2270
14800 1480 4.2699992e-05 0.00023789279 2270
15000 1500 5.3012016e-05 0.00019933234 2270
15200 1520 5.8834197e-05 0.00022407007 2270
15400 1540 5.0899982e-05 0.00029695531 2270
15600 1560 3.0476742e-05 0.00039119066 2270
15800 1580 1.6633264e-05 0.00033770401 2270
16000 1600 1.098906e-05 0.00036684894 2270
16200 1620 1.464848e-05 0.00036449759 2270
16400 1640 1.9598429e-05 0.00021056689 2270
16600 1660 1.2644955e-05 0.00020781781 2270
16800 1680 8.8428553e-06 0.000165 2270
17000 1700 8.8971439e-06 0.00012266475 2270
17200 1720 1.7032781e-05 0.00019873443 2270
17400 1740 1.9448563e-05 0.00025661663 2270
17600 1760 1.3714713e-05 0.000324022 2270
17800 1780 9.1326468e-06 0.00031392513 2270
18000 1800 9.2464802e-06 0.00029729527 2270
18200 1820 1.5553042e-05 0.00027488475 2270
18400 1840 1.4132933e-05 0.00019565459 2270
18600 1860 9.4734832e-06 0.00016716988 2270
18800 1880 5.5115145e-06 0.00013728033 2270
19000 1900 8.268812e-06 0.00015119605 2270
19200 1920 1.2470136e-05 0.00020222131 2270
19400 1940 9.9387775e-06 0.00024503373 2270
19600 1960 5.4241999e-06 0.00026921858 2270
19800 1980 2.7987348e-06 0.00026201267 2270
20000 2000 6.272538e-06 0.00025626323 2270
20200 2020 8.0157781e-06 0.000220139 2270
20400 2040 6.1652093e-06 0.00017089058 2270
20600 2060 2.9967592e-06 0.00014582864 2270
20800 2080 3.016678e-06 0.000148629 2270
21000 2100 7.287645e-06 0.00016486102 2270
21200 2120 8.6905277e-06 0.00020276916 2270
21400 2140 6.8453018e-06 0.00023156153 2270
21600 2160 3.3853799e-06 0.0002432462 2270
21800 2180 4.1241209e-06 0.00022829024 2270
22000 2200 7.0802396e-06 0.00020784823 2270
22200 2220 7.3361691e-06 0.00018114134 2270
22400 2240 5.0764593e-06 0.00014351106 2270
22600 2260 2.7487537e-06 0.00012919872 2270
22800 2280 4.620167e-06 0.00013746218 2270
23000 2300 6.9819357e-06 0.00015985102 2270
23200 2320 6.8923916e-06 0.00018713045 2270
23400 2340 4.1795088e-06 0.00019846682 2270
23600 2360 2.2871028e-06 0.00021068421 2270
23800 2380 3.862046e-06 0.00019553306 2270
24000 2400 5.2448555e-06 0.00017398041 2270
24200 2420 4.7565441e-06 0.00015008142 2270
24400 2440 2.2952135e-06 0.00012747106 2270
24600 2460 2.1575617e-06 0.00012516996 2270
24800 2480 4.1777868e-06 0.0001331902 2270
25000 2500 5.5679133e-06 0.00015504562 2270
25200 2520 4.5758741e-06 0.00017146032 2270
25400 2540 2.3403277e-06 0.00017611666 2270
25600 2560 2.7029302e-06 0.00016850788 2270
25800 2580 4.3601102e-06 0.00015884642 2270
26000 2600 5.2244249e-06 0.00013793898 2270
26200 2620 3.4577672e-06 0.00012395875 2270
26400 2640 2.361577e-06 0.00011600057 2270
26600 2660 2.8515644e-06 0.00011277063 2270
26800 2680 4.0851213e-06 0.0001290832 2270
27000 2700 4.2579644e-06 0.0001476495 2270
27200 2720 2.6593858e-06 0.00015977745 2270
27400 2740 1.990115e-06 0.00015612787 2270
27600 2760 2.6756835e-06 0.00014913772 2270
27800 2780 3.9032806e-06 0.00014014763 2270
28000 2800 3.2729446e-06 0.00012216846 2270
28200 2820 1.9357278e-06 0.00011078621 2270
28400 2840 1.7094832e-06 0.00010910509 2270
28600 2860 2.8731406e-06 0.00011179644 2270
28800 2880 3.7062354e-06 0.00012254091 2270
29000 2900 2.7844262e-06 0.00013060331 2270
29200 2920 1.7680655e-06 0.00013797514 2270
29400 2940 1.706873e-06 0.0001350685 2270
29600 2960 2.8764562e-06 0.00012428508 2270
29800 2980 3.1502029e-06 0.00011456718 2270
30000 3000 2.1833409e-06 0.00010317469 2270
Loop time of 165.611 on 4 procs for 30000 steps with 2270 atoms
Performance: 1565111.240 tau/day, 181.147 timesteps/s, 411.204 katom-step/s
99.7% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.63183 | 21.226 | 42.266 | 444.6 | 12.82
Bond | 0.095073 | 0.17799 | 0.27877 | 17.0 | 0.11
Neigh | 2.0745 | 2.0781 | 2.0822 | 0.2 | 1.25
Comm | 0.32024 | 0.38703 | 0.45564 | 8.1 | 0.23
Output | 0.60459 | 0.76798 | 0.93724 | 18.6 | 0.46
Modify | 119.85 | 140.76 | 161.36 | 172.2 | 85.00
Other | | 0.2124 | | | 0.13
Nlocal: 567.5 ave 1139 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 75.5 ave 152 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 9238.25 ave 18490 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
FullNghs: 17945 ave 35917 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 71780
Ave neighs/atom = 31.621145
Ave special neighs/atom = 0.22026432
Neighbor list builds = 2071
Dangerous builds = 0
Total wall time: 0:02:45

658
examples/rheo/ice-cubes/square.mol Normal file
View File

@ -0,0 +1,658 @@
#Made with create_mol.py
100 atoms
342 bonds
Coords
#ID x y z
1 -4 -4 0
2 -3 -4 0
3 -2 -4 0
4 -1 -4 0
5 0 -4 0
6 1 -4 0
7 2 -4 0
8 3 -4 0
9 4 -4 0
10 5 -4 0
11 -4 -3 0
12 -3 -3 0
13 -2 -3 0
14 -1 -3 0
15 0 -3 0
16 1 -3 0
17 2 -3 0
18 3 -3 0
19 4 -3 0
20 5 -3 0
21 -4 -2 0
22 -3 -2 0
23 -2 -2 0
24 -1 -2 0
25 0 -2 0
26 1 -2 0
27 2 -2 0
28 3 -2 0
29 4 -2 0
30 5 -2 0
31 -4 -1 0
32 -3 -1 0
33 -2 -1 0
34 -1 -1 0
35 0 -1 0
36 1 -1 0
37 2 -1 0
38 3 -1 0
39 4 -1 0
40 5 -1 0
41 -4 0 0
42 -3 0 0
43 -2 0 0
44 -1 0 0
45 0 0 0
46 1 0 0
47 2 0 0
48 3 0 0
49 4 0 0
50 5 0 0
51 -4 1 0
52 -3 1 0
53 -2 1 0
54 -1 1 0
55 0 1 0
56 1 1 0
57 2 1 0
58 3 1 0
59 4 1 0
60 5 1 0
61 -4 2 0
62 -3 2 0
63 -2 2 0
64 -1 2 0
65 0 2 0
66 1 2 0
67 2 2 0
68 3 2 0
69 4 2 0
70 5 2 0
71 -4 3 0
72 -3 3 0
73 -2 3 0
74 -1 3 0
75 0 3 0
76 1 3 0
77 2 3 0
78 3 3 0
79 4 3 0
80 5 3 0
81 -4 4 0
82 -3 4 0
83 -2 4 0
84 -1 4 0
85 0 4 0
86 1 4 0
87 2 4 0
88 3 4 0
89 4 4 0
90 5 4 0
91 -4 5 0
92 -3 5 0
93 -2 5 0
94 -1 5 0
95 0 5 0
96 1 5 0
97 2 5 0
98 3 5 0
99 4 5 0
100 5 5 0
Types
#ID type
1 1
2 1
3 1
4 1
5 1
6 1
7 1
8 1
9 1
10 1
11 1
12 1
13 1
14 1
15 1
16 1
17 1
18 1
19 1
20 1
21 1
22 1
23 1
24 1
25 1
26 1
27 1
28 1
29 1
30 1
31 1
32 1
33 1
34 1
35 1
36 1
37 1
38 1
39 1
40 1
41 1
42 1
43 1
44 1
45 1
46 1
47 1
48 1
49 1
50 1
51 1
52 1
53 1
54 1
55 1
56 1
57 1
58 1
59 1
60 1
61 1
62 1
63 1
64 1
65 1
66 1
67 1
68 1
69 1
70 1
71 1
72 1
73 1
74 1
75 1
76 1
77 1
78 1
79 1
80 1
81 1
82 1
83 1
84 1
85 1
86 1
87 1
88 1
89 1
90 1
91 1
92 1
93 1
94 1
95 1
96 1
97 1
98 1
99 1
100 1
Masses
#ID mass
1 1
2 1
3 1
4 1
5 1
6 1
7 1
8 1
9 1
10 1
11 1
12 1
13 1
14 1
15 1
16 1
17 1
18 1
19 1
20 1
21 1
22 1
23 1
24 1
25 1
26 1
27 1
28 1
29 1
30 1
31 1
32 1
33 1
34 1
35 1
36 1
37 1
38 1
39 1
40 1
41 1
42 1
43 1
44 1
45 1
46 1
47 1
48 1
49 1
50 1
51 1
52 1
53 1
54 1
55 1
56 1
57 1
58 1
59 1
60 1
61 1
62 1
63 1
64 1
65 1
66 1
67 1
68 1
69 1
70 1
71 1
72 1
73 1
74 1
75 1
76 1
77 1
78 1
79 1
80 1
81 1
82 1
83 1
84 1
85 1
86 1
87 1
88 1
89 1
90 1
91 1
92 1
93 1
94 1
95 1
96 1
97 1
98 1
99 1
100 1
Bonds
#ID type atom1 atom2
1 1 1 2
2 1 1 11
3 1 1 12
4 1 2 3
5 1 2 11
6 1 2 12
7 1 2 13
8 1 3 4
9 1 3 12
10 1 3 13
11 1 3 14
12 1 4 5
13 1 4 13
14 1 4 14
15 1 4 15
16 1 5 6
17 1 5 14
18 1 5 15
19 1 5 16
20 1 6 7
21 1 6 15
22 1 6 16
23 1 6 17
24 1 7 8
25 1 7 16
26 1 7 17
27 1 7 18
28 1 8 9
29 1 8 17
30 1 8 18
31 1 8 19
32 1 9 10
33 1 9 18
34 1 9 19
35 1 9 20
36 1 10 19
37 1 10 20
38 1 11 21
39 1 11 12
40 1 11 22
41 1 12 21
42 1 12 13
43 1 12 22
44 1 12 23
45 1 13 22
46 1 13 23
47 1 13 14
48 1 13 24
49 1 14 23
50 1 14 24
51 1 14 15
52 1 14 25
53 1 15 24
54 1 15 16
55 1 15 25
56 1 15 26
57 1 16 25
58 1 16 26
59 1 16 17
60 1 16 27
61 1 17 26
62 1 17 18
63 1 17 27
64 1 17 28
65 1 18 27
66 1 18 28
67 1 18 19
68 1 18 29
69 1 19 28
70 1 19 29
71 1 19 20
72 1 19 30
73 1 20 29
74 1 20 30
75 1 21 22
76 1 21 31
77 1 21 32
78 1 22 23
79 1 22 31
80 1 22 32
81 1 22 33
82 1 23 24
83 1 23 32
84 1 23 33
85 1 23 34
86 1 24 25
87 1 24 33
88 1 24 34
89 1 24 35
90 1 25 26
91 1 25 34
92 1 25 35
93 1 25 36
94 1 26 27
95 1 26 35
96 1 26 36
97 1 26 37
98 1 27 28
99 1 27 36
100 1 27 37
101 1 27 38
102 1 28 29
103 1 28 37
104 1 28 38
105 1 28 39
106 1 29 30
107 1 29 38
108 1 29 39
109 1 29 40
110 1 30 39
111 1 30 40
112 1 31 32
113 1 31 41
114 1 31 42
115 1 32 33
116 1 32 41
117 1 32 42
118 1 32 43
119 1 33 34
120 1 33 42
121 1 33 43
122 1 33 44
123 1 34 35
124 1 34 43
125 1 34 44
126 1 34 45
127 1 35 36
128 1 35 44
129 1 35 45
130 1 35 46
131 1 36 37
132 1 36 45
133 1 36 46
134 1 36 47
135 1 37 38
136 1 37 46
137 1 37 47
138 1 37 48
139 1 38 39
140 1 38 47
141 1 38 48
142 1 38 49
143 1 39 40
144 1 39 48
145 1 39 49
146 1 39 50
147 1 40 49
148 1 40 50
149 1 41 51
150 1 41 42
151 1 41 52
152 1 42 51
153 1 42 43
154 1 42 52
155 1 42 53
156 1 43 52
157 1 43 53
158 1 43 44
159 1 43 54
160 1 44 53
161 1 44 54
162 1 44 45
163 1 44 55
164 1 45 54
165 1 45 46
166 1 45 55
167 1 45 56
168 1 46 55
169 1 46 56
170 1 46 47
171 1 46 57
172 1 47 56
173 1 47 48
174 1 47 57
175 1 47 58
176 1 48 57
177 1 48 58
178 1 48 49
179 1 48 59
180 1 49 58
181 1 49 59
182 1 49 50
183 1 49 60
184 1 50 59
185 1 50 60
186 1 51 52
187 1 51 61
188 1 51 62
189 1 52 53
190 1 52 61
191 1 52 62
192 1 52 63
193 1 53 54
194 1 53 62
195 1 53 63
196 1 53 64
197 1 54 55
198 1 54 63
199 1 54 64
200 1 54 65
201 1 55 56
202 1 55 64
203 1 55 65
204 1 55 66
205 1 56 57
206 1 56 65
207 1 56 66
208 1 56 67
209 1 57 58
210 1 57 66
211 1 57 67
212 1 57 68
213 1 58 59
214 1 58 67
215 1 58 68
216 1 58 69
217 1 59 60
218 1 59 68
219 1 59 69
220 1 59 70
221 1 60 69
222 1 60 70
223 1 61 71
224 1 61 62
225 1 61 72
226 1 62 71
227 1 62 63
228 1 62 72
229 1 62 73
230 1 63 72
231 1 63 73
232 1 63 64
233 1 63 74
234 1 64 73
235 1 64 74
236 1 64 65
237 1 64 75
238 1 65 74
239 1 65 66
240 1 65 75
241 1 65 76
242 1 66 75
243 1 66 76
244 1 66 67
245 1 66 77
246 1 67 76
247 1 67 68
248 1 67 77
249 1 67 78
250 1 68 77
251 1 68 78
252 1 68 69
253 1 68 79
254 1 69 78
255 1 69 79
256 1 69 70
257 1 69 80
258 1 70 79
259 1 70 80
260 1 71 72
261 1 71 81
262 1 71 82
263 1 72 73
264 1 72 81
265 1 72 82
266 1 72 83
267 1 73 74
268 1 73 82
269 1 73 83
270 1 73 84
271 1 74 75
272 1 74 83
273 1 74 84
274 1 74 85
275 1 75 76
276 1 75 84
277 1 75 85
278 1 75 86
279 1 76 77
280 1 76 85
281 1 76 86
282 1 76 87
283 1 77 78
284 1 77 86
285 1 77 87
286 1 77 88
287 1 78 79
288 1 78 87
289 1 78 88
290 1 78 89
291 1 79 80
292 1 79 88
293 1 79 89
294 1 79 90
295 1 80 89
296 1 80 90
297 1 81 82
298 1 81 91
299 1 81 92
300 1 82 83
301 1 82 91
302 1 82 92
303 1 82 93
304 1 83 84
305 1 83 92
306 1 83 93
307 1 83 94
308 1 84 85
309 1 84 93
310 1 84 94
311 1 84 95
312 1 85 86
313 1 85 94
314 1 85 95
315 1 85 96
316 1 86 87
317 1 86 95
318 1 86 96
319 1 86 97
320 1 87 88
321 1 87 96
322 1 87 97
323 1 87 98
324 1 88 89
325 1 88 97
326 1 88 98
327 1 88 99
328 1 89 90
329 1 89 98
330 1 89 99
331 1 89 100
332 1 90 99
333 1 90 100
334 1 91 92
335 1 92 93
336 1 93 94
337 1 94 95
338 1 95 96
339 1 96 97
340 1 97 98
341 1 98 99
342 1 99 100

102
examples/rheo/oxidation/in.rheo.oxidation Normal file
View File

@ -0,0 +1,102 @@
# ------ 2D oxidizing bar ------ #
dimension 2
units lj
atom_style hybrid rheo/thermal bond
boundary m m p
comm_modify vel yes
newton off
region box block -60 60 0 80 -0.01 0.01 units box
create_box 3 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50
region lbar block -15 0 3 80 EDGE EDGE units box
region rbar block 0 15 3 80 EDGE EDGE units box
region bar union 2 lbar rbar
region floor block EDGE EDGE EDGE 3.0 EDGE EDGE units box
lattice hex 1.0
create_atoms 1 region bar
create_atoms 3 region floor
set region rbar type 2
group bar type 1 2
group rbar type 2
group floor type 3
set group all sph/e 0.0
set group all rheo/status 1
# ------ Model parameters ------#
variable cut equal 3.0
variable n equal 1.0
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable zeta equal 0.05
variable kappa equal 0.1*${rho0}/${mp}
variable dt_max equal 0.1*${cut}/${cs}/3
variable eta equal 0.05
variable Cv equal 1.0
variable L equal 0.1
variable Tf equal 1.0
mass * ${mp}
timestep 0.1
pair_style hybrid/overlay rheo ${cut} artificial/visc ${zeta} rheo/solid
pair_coeff * * rheo
pair_coeff * * rheo/solid 1.0 1.0 1.0
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
create_bonds many bar bar 1 0 1.5
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
bond_style hybrid bpm/spring rheo/shell t/form 100
bond_coeff 1 bpm/spring 1.0 1.0 1.0
bond_coeff 2 rheo/shell 0.2 0.2 0.1
# ------ Apply dynamics ------#
# Note: surface detection is not performed on solid bodies, so cannot use surface property
compute coord all rheo/property/atom coordination
variable surf atom c_coord<22
group surf dynamic all var surf every 10
fix 1 all rheo ${cut} quintic 0 &
thermal &
shift &
surface/detection coordination 22 8
fix 2 all rheo/viscosity * constant ${eta}
fix 3 all rheo/pressure * linear
fix 4 all rheo/thermal conductivity * constant ${kappa} &
specific/heat * constant ${Cv} &
Tfreeze * constant ${Tf} &
latent/heat * constant ${L} &
react 1.5 1
fix 5 rbar rheo/oxidation 1.5 2 1.0
fix 6 all wall/harmonic ylo EDGE 2.0 1.0 1.0
fix 7 all gravity 5e-5 vector 0 -1 0
fix 8 floor setforce 0.0 0.0 0.0
fix 9 surf add/heat linear 1.1 0.05
fix 10 floor add/heat constant 0 overwrite yes # fix the temperature of the floor
fix 11 all enforce2d
compute surf all rheo/property/atom surface
compute rho all rheo/property/atom rho
compute phase all rheo/property/atom phase
compute temp all rheo/property/atom temperature
compute eng all rheo/property/atom energy
compute nbond_shell all rheo/property/atom nbond/shell
compute nbond_solid all nbond/atom bond/type 1
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press atoms
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_phase c_temp c_eng c_nbond_solid c_nbond_shell c_rho c_surf
run 40000

488
examples/rheo/oxidation/log.17Apr2024.oxidation.g++.4 Normal file
View File

@ -0,0 +1,488 @@
LAMMPS (17 Apr 2024 - Development - patch_5May2020-18508-g3c0eaf6870-modified)
# ------ 2D oxidizing bar ------ #
dimension 2
units lj
atom_style hybrid rheo/thermal bond
boundary m m p
comm_modify vel yes
newton off
region box block -60 60 0 80 -0.01 0.01 units box
create_box 3 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50
Created orthogonal box = (-60 0 -0.01) to (60 80 0.01)
2 by 2 by 1 MPI processor grid
region lbar block -15 0 3 80 EDGE EDGE units box
region rbar block 0 15 3 80 EDGE EDGE units box
region bar union 2 lbar rbar
region floor block EDGE EDGE EDGE 3.0 EDGE EDGE units box
lattice hex 1.0
Lattice spacing in x,y,z = 1.0745699 1.8612097 1.0745699
create_atoms 1 region bar
Created 2255 atoms
using lattice units in orthogonal box = (-60 0 -0.01) to (60 80 0.01)
create_atoms CPU = 0.001 seconds
create_atoms 3 region floor
Created 446 atoms
using lattice units in orthogonal box = (-60 0 -0.01) to (60 80 0.01)
create_atoms CPU = 0.000 seconds
set region rbar type 2
Setting atom values ...
1148 settings made for type
group bar type 1 2
2255 atoms in group bar
group rbar type 2
1148 atoms in group rbar
group floor type 3
446 atoms in group floor
set group all sph/e 0.0
Setting atom values ...
2701 settings made for sph/e
set group all rheo/status 1
Setting atom values ...
2701 settings made for rheo/status
# ------ Model parameters ------#
variable cut equal 3.0
variable n equal 1.0
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable mp equal 1/${n}
variable mp equal 1/1
variable zeta equal 0.05
variable kappa equal 0.1*${rho0}/${mp}
variable kappa equal 0.1*1/${mp}
variable kappa equal 0.1*1/1
variable dt_max equal 0.1*${cut}/${cs}/3
variable dt_max equal 0.1*3/${cs}/3
variable dt_max equal 0.1*3/1/3
variable eta equal 0.05
variable Cv equal 1.0
variable L equal 0.1
variable Tf equal 1.0
mass * ${mp}
mass * 1
timestep 0.1
pair_style hybrid/overlay rheo ${cut} artificial/visc ${zeta} rheo/solid
pair_style hybrid/overlay rheo 3 artificial/visc ${zeta} rheo/solid
pair_style hybrid/overlay rheo 3 artificial/visc 0.05 rheo/solid
pair_coeff * * rheo
pair_coeff * * rheo/solid 1.0 1.0 1.0
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
0 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.000 seconds
create_bonds many bar bar 1 0 1.5
Generated 0 of 3 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 = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 73 49 1
3 neighbor lists, perpetual/occasional/extra = 2 1 0
(1) command create_bonds, occasional
attributes: full, newton off
pair build: full/bin
stencil: full/bin/2d
bin: standard
(2) pair rheo, perpetual
attributes: half, newton off
pair build: half/bin/newtoff
stencil: full/bin/2d
bin: standard
(3) pair rheo/solid, perpetual, trim from (2)
attributes: half, newton off, cut 1.3
pair build: trim
stencil: none
bin: none
Added 6547 bonds, new total = 6547
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 1 1
special bond factors coul: 0 1 1
6 = max # of 1-2 neighbors
101 = max # of special neighbors
special bonds CPU = 0.000 seconds
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
bond_style hybrid bpm/spring rheo/shell t/form 100
bond_coeff 1 bpm/spring 1.0 1.0 1.0
bond_coeff 2 rheo/shell 0.2 0.2 0.1
# ------ Apply dynamics ------#
# Note: surface detection is not performed on solid bodies, so cannot use surface property
compute coord all rheo/property/atom coordination
variable surf atom c_coord<22
group surf dynamic all var surf every 10
dynamic group surf defined
fix 1 all rheo ${cut} quintic 0 thermal shift surface/detection coordination 22 8
fix 1 all rheo 3 quintic 0 thermal shift surface/detection coordination 22 8
fix 2 all rheo/viscosity * constant ${eta}
fix 2 all rheo/viscosity * constant 0.05
fix 3 all rheo/pressure * linear
fix 4 all rheo/thermal conductivity * constant ${kappa} specific/heat * constant ${Cv} Tfreeze * constant ${Tf} latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.1 specific/heat * constant ${Cv} Tfreeze * constant ${Tf} latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.1 specific/heat * constant 1 Tfreeze * constant ${Tf} latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.1 specific/heat * constant 1 Tfreeze * constant 1 latent/heat * constant ${L} react 1.5 1
fix 4 all rheo/thermal conductivity * constant 0.1 specific/heat * constant 1 Tfreeze * constant 1 latent/heat * constant 0.1 react 1.5 1
fix 5 rbar rheo/oxidation 1.5 2 1.0
fix 6 all wall/harmonic ylo EDGE 2.0 1.0 1.0
fix 7 all gravity 5e-5 vector 0 -1 0
fix 8 floor setforce 0.0 0.0 0.0
fix 9 surf add/heat linear 1.1 0.05
fix 10 floor add/heat constant 0 overwrite yes # fix the temperature of the floor
fix 11 all enforce2d
compute surf all rheo/property/atom surface
compute rho all rheo/property/atom rho
compute phase all rheo/property/atom phase
compute status all rheo/property/atom status
compute temp all rheo/property/atom temperature
compute eng all rheo/property/atom energy
compute nbond_shell all rheo/property/atom nbond/shell
compute nbond_solid all nbond/atom bond/type 1
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press atoms
dump 1 all custom 200 atomDump id type x y vx vy fx fy c_phase c_temp c_eng c_nbond_solid c_nbond_shell c_rho c_surf c_status
run 40000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- BPM bond style: doi:10.1039/D3SM01373A
@Article{Clemmer2024,
author = {Clemmer, Joel T. and Monti, Joseph M. and Lechman, Jeremy B.},
title = {A soft departure from jamming: the compaction of deformable
granular matter under high pressures},
journal = {Soft Matter},
year = 2024,
volume = 20,
number = 8,
pages = {1702--1718}
}
- @article{PalermoInPrep,
journal = {in prep},
title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},
year = {2024},
author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},
}
- @article{ApplMathModel.130.310,
title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},
journal = {Applied Mathematical Modelling},
volume = {130},
pages = {310-326},
year = {2024},
issn = {0307-904X},
doi = {https://doi.org/10.1016/j.apm.2024.02.027},
author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 3 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 = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 73 49 1
8 neighbor lists, perpetual/occasional/extra = 7 1 0
(1) pair rheo, perpetual, half/full from (3)
attributes: half, newton off
pair build: halffull/newtoff
stencil: none
bin: none
(2) pair rheo/solid, perpetual, trim from (4)
attributes: half, newton off, cut 1.3
pair build: trim
stencil: none
bin: none
(3) compute RHEO/KERNEL, perpetual
attributes: full, newton off
pair build: full/bin
stencil: full/bin/2d
bin: standard
(4) compute RHEO/GRAD, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(5) compute RHEO/VSHIFT, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(6) compute RHEO/SURFACE, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(7) fix rheo/thermal, occasional, trim from (4)
attributes: half, newton off, cut 3
pair build: trim
stencil: none
bin: none
(8) fix rheo/oxidation, perpetual, trim from (3)
attributes: full, newton off, cut 1.8
pair build: trim
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 25.96 | 25.96 | 25.96 Mbytes
Step Time KinEng Press Atoms
0 0 0 0 2701
200 20 4.1743799e-07 1.1743617e-07 2701
400 40 1.6697519e-06 4.6974469e-07 2701
600 60 3.7127333e-06 1.0646825e-05 2701
800 80 4.6683656e-06 0.00015182605 2701
1000 100 4.7368707e-06 0.00028128761 2701
1200 120 3.4384322e-06 0.00045913378 2701
1400 140 1.4119866e-06 0.00055627091 2701
1600 160 4.4114517e-07 0.00058247308 2701
1800 180 4.8289229e-07 0.0005510948 2701
2000 200 1.8494183e-06 0.00048386222 2701
2200 220 3.3319816e-06 0.00037903264 2701
2400 240 3.8128922e-06 0.00024115906 2701
2600 260 3.1943401e-06 9.727407e-05 2701
2800 280 1.6172816e-06 -2.632162e-05 2701
3000 300 3.6100709e-07 -8.5761867e-05 2701
3200 320 1.4745502e-07 -5.9204127e-05 2701
3400 340 8.3369782e-07 8.8312464e-07 2701
3600 360 2.0484052e-06 5.8521477e-05 2701
3800 380 3.1639387e-06 0.0001685663 2701
4000 400 3.1692907e-06 0.00026875988 2701
4200 420 2.391933e-06 0.00038621787 2701
4400 440 1.1964404e-06 0.00048901286 2701
4600 460 4.0508824e-07 0.00051863639 2701
4800 480 5.4908507e-07 0.00049263754 2701
5000 500 1.3139665e-06 0.00041984264 2701
5200 520 2.1939161e-06 0.00033095351 2701
5400 540 2.3687031e-06 0.00022422981 2701
5600 560 1.8280882e-06 0.00011544328 2701
5800 580 8.8610517e-07 2.9307791e-05 2701
6000 600 2.0989359e-07 -1.7340941e-05 2701
6200 620 2.8658301e-07 -8.1237835e-06 2701
6400 640 9.7636239e-07 4.3755922e-05 2701
6600 660 1.891303e-06 0.0001185719 2701
6800 680 2.4149904e-06 0.00020830273 2701
7000 700 2.3174953e-06 0.00030114767 2701
7200 720 1.7918612e-06 0.00037821537 2701
7400 740 1.2114987e-06 0.0004233475 2701
7600 760 9.9661553e-07 0.00042958263 2701
7800 780 1.1552559e-06 0.00039944618 2701
8000 800 1.5249138e-06 0.00034034478 2701
8200 820 1.7453861e-06 0.00026826463 2701
8400 840 1.6259021e-06 0.00019131768 2701
8600 860 1.2612805e-06 0.0001162957 2701
8800 880 8.6964518e-07 7.1771506e-05 2701
9000 900 7.6892472e-07 5.6170687e-05 2701
9200 920 1.0780045e-06 7.1925995e-05 2701
9400 940 1.6514902e-06 0.00011635293 2701
9600 960 2.1891377e-06 0.00017599885 2701
9800 980 2.4551701e-06 0.00024127934 2701
10000 1000 2.4277051e-06 0.00029918622 2701
10200 1020 2.2655987e-06 0.00034067996 2701
10400 1040 2.1767207e-06 0.00035598133 2701
10600 1060 2.2796719e-06 0.00034359076 2701
10800 1080 2.4884225e-06 0.00030749714 2701
11000 1100 2.6387215e-06 0.00025725198 2701
11200 1120 2.5968908e-06 0.00020170699 2701
11400 1140 2.4108931e-06 0.00015185858 2701
11600 1160 2.2375166e-06 0.00011800349 2701
11800 1180 2.2407196e-06 0.00010646971 2701
12000 1200 2.4845263e-06 0.00011817498 2701
12200 1220 2.8733204e-06 0.00015013186 2701
12400 1240 3.2437087e-06 0.00019211975 2701
12600 1260 3.4732728e-06 0.00023620276 2701
12800 1280 3.5836611e-06 0.00027352269 2701
13000 1300 3.6592211e-06 0.00029533734 2701
13200 1320 3.782506e-06 0.00030032559 2701
13400 1340 3.9807086e-06 0.00028395722 2701
13600 1360 4.2023176e-06 0.00025390325 2701
13800 1380 4.3559781e-06 0.00021794236 2701
14000 1400 4.4273371e-06 0.00018026034 2701
14200 1420 4.49867e-06 0.0001526569 2701
14400 1440 4.6591574e-06 0.00013707051 2701
14600 1460 4.9589583e-06 0.00013803875 2701
14800 1480 5.3859375e-06 0.00015455425 2701
15000 1500 5.8639557e-06 0.00017954785 2701
15200 1520 6.3075561e-06 0.0002084257 2701
15400 1540 6.7022179e-06 0.0002347669 2701
15600 1560 7.0789688e-06 0.00025020766 2701
15800 1580 7.4734777e-06 0.00025394845 2701
16000 1600 7.8884743e-06 0.00024571725 2701
16200 1620 8.3224059e-06 0.00022706648 2701
16400 1640 8.7337783e-06 0.00020320706 2701
16600 1660 9.1454649e-06 0.00017824346 2701
16800 1680 9.5948793e-06 0.00015961835 2701
17000 1700 1.0106407e-05 0.00015135471 2701
17200 1720 1.0707273e-05 0.00015166884 2701
17400 1740 1.1392597e-05 0.0001645916 2701
17600 1760 1.2118829e-05 0.00018119729 2701
17800 1780 1.2846056e-05 0.0002003616 2701
18000 1800 1.3555288e-05 0.00021585952 2701
18200 1820 1.4301024e-05 0.00022290158 2701
18400 1840 1.5089217e-05 0.00021970192 2701
18600 1860 1.5902351e-05 0.00020911128 2701
18800 1880 1.6753175e-05 0.00019278718 2701
19000 1900 1.7602996e-05 0.00017584076 2701
19200 1920 1.8479378e-05 0.00016206226 2701
19400 1940 1.9421603e-05 0.00015575677 2701
19600 1960 2.0477421e-05 0.00015687558 2701
19800 1980 2.1617288e-05 0.00016424998 2701
20000 2000 2.2814347e-05 0.00017466664 2701
20200 2020 2.4029097e-05 0.00018647149 2701
20400 2040 2.5255953e-05 0.00019516077 2701
20600 2060 2.649418e-05 0.00019906384 2701
20800 2080 2.7755897e-05 0.00019630586 2701
21000 2100 2.9067854e-05 0.00018674721 2701
21200 2120 3.0396477e-05 0.0001758048 2701
21400 2140 3.1759719e-05 0.00016782801 2701
21600 2160 3.3193597e-05 0.00016324138 2701
21800 2180 3.4729384e-05 0.00016124274 2701
22000 2200 3.6367594e-05 0.00016437457 2701
22200 2220 3.8095131e-05 0.00017015573 2701
22400 2240 3.9867003e-05 0.00017649465 2701
22600 2260 4.169511e-05 0.00018111374 2701
22800 2280 4.3566134e-05 0.00018104136 2701
23000 2300 4.5461538e-05 0.00017822707 2701
23200 2320 4.7377333e-05 0.00017285066 2701
23400 2340 4.9354403e-05 0.00016826524 2701
23600 2360 5.1399791e-05 0.00016517913 2701
23800 2380 5.3510931e-05 0.00016299649 2701
24000 2400 5.5681048e-05 0.00016256674 2701
24200 2420 5.7902429e-05 0.00016513449 2701
24400 2440 6.0216049e-05 0.00016895109 2701
24600 2460 6.270982e-05 0.00016946227 2701
24800 2480 6.5390117e-05 0.00016589426 2701
25000 2500 6.8121899e-05 0.00016241676 2701
25200 2520 7.0947331e-05 0.00015624292 2701
25400 2540 7.4304148e-05 0.0001449537 2701
25600 2560 7.7745077e-05 0.00013179658 2701
25800 2580 8.0739829e-05 0.00013098838 2701
26000 2600 8.3827874e-05 0.00014278841 2701
26200 2620 8.7060677e-05 0.00015381649 2701
26400 2640 9.0266508e-05 0.00016130999 2701
26600 2660 9.3339049e-05 0.00016908268 2701
26800 2680 9.6347013e-05 0.00016771087 2701
27000 2700 9.9294711e-05 0.00016577315 2701
27200 2720 0.00010230007 0.0001670893 2701
27400 2740 0.00010547172 0.00016569077 2701
27600 2760 0.00010872426 0.00016506303 2701
27800 2780 0.00011201844 0.00016482702 2701
28000 2800 0.00011532129 0.00016694886 2701
28200 2820 0.00011869854 0.00016163005 2701
28400 2840 0.00012209747 0.00015339281 2701
28600 2860 0.00012549322 0.00014765883 2701
28800 2880 0.00012898685 0.00014241765 2701
29000 2900 0.00013259039 0.00014215724 2701
29200 2920 0.00013628209 0.00014881155 2701
29400 2940 0.00014001213 0.00015671333 2701
29600 2960 0.00014379216 0.00016446215 2701
29800 2980 0.00014764687 0.0001639602 2701
30000 3000 0.00015142301 0.00015664816 2701
30200 3020 0.00015496407 0.00015545099 2701
30400 3040 0.00015797338 0.00015368625 2701
30600 3060 0.00016042141 0.00015679918 2701
30800 3080 0.00016244716 0.00016093678 2701
31000 3100 0.00016202247 0.00016066954 2701
31200 3120 0.0001613312 0.00015932059 2701
31400 3140 0.00016274961 0.00015988567 2701
31600 3160 0.00016541518 0.00015724809 2701
31800 3180 0.00016809362 0.00015498827 2701
32000 3200 0.00017067801 0.00014830489 2701
32200 3220 0.00017333906 0.00014371345 2701
32400 3240 0.0001759011 0.00014421259 2701
32600 3260 0.00017849952 0.00014228443 2701
32800 3280 0.00017801812 0.00014117391 2701
33000 3300 0.00017718857 0.00014644675 2701
33200 3320 0.00017833666 0.0001291286 2701
33400 3340 0.000178576 0.00014878558 2701
33600 3360 0.00017846711 0.00013905481 2701
33800 3380 0.00017822937 0.00015535996 2701
34000 3400 0.00017899663 0.00016094303 2701
34200 3420 0.00017924661 0.00015017553 2701
34400 3440 0.00018024855 0.00014723549 2701
34600 3460 0.00018143865 0.00013903131 2701
34800 3480 0.00018258173 0.00013722112 2701
35000 3500 0.00018404873 0.00014675949 2701
35200 3520 0.00018538521 0.00015108242 2701
35400 3540 0.00018669649 0.00014564852 2701
35600 3560 0.00018814608 0.00013762161 2701
35800 3580 0.00018967415 0.00014602307 2701
36000 3600 0.00019146735 0.000126909 2701
36200 3620 0.00019414036 0.00012384379 2701
36400 3640 0.00019613057 0.00011059573 2701
36600 3660 0.00019897104 0.00013621801 2701
36800 3680 0.00020169688 0.00013665462 2701
37000 3700 0.00020447655 0.00013929258 2701
37200 3720 0.00020711105 0.0001363895 2701
37400 3740 0.00021077854 0.00013610672 2701
37600 3760 0.00021303084 0.00015051235 2701
37800 3780 0.00021619561 0.00012664801 2701
38000 3800 0.0002194018 0.00012808247 2701
38200 3820 0.00022242646 0.0001360174 2701
38400 3840 0.00022531568 0.00013311221 2701
38600 3860 0.00022821731 0.00013523939 2701
38800 3880 0.000231228 0.00014090695 2701
39000 3900 0.00023404038 0.00013661835 2701
39200 3920 0.00023755044 0.00013659469 2701
39400 3940 0.00024009059 0.00012097907 2701
39600 3960 0.0002432098 9.7877876e-05 2701
39800 3980 0.00024475294 0.0001164688 2701
40000 4000 0.00024171274 0.00012432219 2701
Loop time of 192.659 on 4 procs for 40000 steps with 2701 atoms
Performance: 1793840.118 tau/day, 207.620 timesteps/s, 560.783 katom-step/s
99.6% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 16.881 | 24.402 | 30.74 | 114.6 | 12.67
Bond | 1.1126 | 1.8917 | 2.6935 | 43.3 | 0.98
Neigh | 35.387 | 35.508 | 35.625 | 1.5 | 18.43
Comm | 1.5499 | 1.6694 | 1.8006 | 7.4 | 0.87
Output | 0.99755 | 1.0072 | 1.0165 | 0.8 | 0.52
Modify | 120.6 | 127.43 | 135.54 | 54.8 | 66.14
Other | | 0.7553 | | | 0.39
Nlocal: 675.25 ave 1373 max 7 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 103 ave 163 max 50 min
Histogram: 2 0 0 0 0 0 0 0 1 1
Neighs: 10509 ave 21592 max 126 min
Histogram: 2 0 0 0 0 0 0 0 0 2
FullNghs: 20367 ave 41981 max 141 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 81468
Ave neighs/atom = 30.162162
Ave special neighs/atom = 1.6593854
Neighbor list builds = 39932
Dangerous builds = 0
Total wall time: 0:03:12

75
examples/rheo/poiseuille/in.rheo.poiseuille Normal file
View File

@ -0,0 +1,75 @@
# ------ 2D Poiseuille flow ------ #
dimension 2
units lj
atom_style rheo
boundary p p p
comm_modify vel yes
# ------ Create simulation box ------ #
variable n equal 1.0
variable cut equal 3.0
region box block 0 20 -10 10 -0.01 0.01
create_box 2 box
lattice sq ${n}
region inner block INF INF -7.5 7.5 INF INF units box
region walls block INF INF -7.5 7.5 INF INF units box side out
create_atoms 2 region walls
create_atoms 1 region inner
group fluid type 1
group rig type 2
displace_atoms fluid random 0.1 0.1 0 135414 units box
# ------ Model parameters ------ #
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable zeta equal 1.0
variable kappa equal 1.0*${rho0}/${mp}
variable fext equal 1e-4/${n}
variable dt_max equal 0.1*${cut}/${cs}/3
variable Dr equal 0.05*${cut}*${cs}
variable eta equal 0.1
variable gd0 equal 5e-4
variable npow equal 0.5
variable K equal 0.001
mass * ${mp}
set group all rheo/rho ${rho0}
set group all rheo/status 0
set group rig rheo/status 1
timestep ${dt_max}
pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
pair_coeff * *
# ------ Fixes & computes ------ #
fix 1 all rheo ${cut} quintic 0 shift
fix 2 all rheo/viscosity * constant ${eta}
#fix 2 all rheo/viscosity * power ${eta} ${gd0} ${K} ${npow}
fix 3 all rheo/pressure * linear
fix 4 rig setforce 0.0 0.0 0.0
fix 5 fluid addforce ${fext} 0.0 0.0
fix 6 all enforce2d
compute rho all rheo/property/atom rho
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_rho
run 20000

288
examples/rheo/poiseuille/log.17Apr2024.poiseuille.g++.4 Normal file
View File

@ -0,0 +1,288 @@
LAMMPS (17 Apr 2024 - Development - patch_5May2020-18508-g3c0eaf6870-modified)
# ------ 2D Poiseuille flow ------ #
dimension 2
units lj
atom_style rheo
boundary p p p
comm_modify vel yes
# ------ Create simulation box ------ #
variable n equal 1.0
variable cut equal 3.0
region box block 0 20 -10 10 -0.01 0.01
create_box 2 box
Created orthogonal box = (0 -10 -0.01) to (20 10 0.01)
2 by 2 by 1 MPI processor grid
lattice sq ${n}
lattice sq 1
Lattice spacing in x,y,z = 1 1 1
region inner block INF INF -7.5 7.5 INF INF units box
region walls block INF INF -7.5 7.5 INF INF units box side out
create_atoms 2 region walls
Created 100 atoms
using lattice units in orthogonal box = (0 -10 -0.01) to (20 10 0.01)
create_atoms CPU = 0.000 seconds
create_atoms 1 region inner
Created 300 atoms
using lattice units in orthogonal box = (0 -10 -0.01) to (20 10 0.01)
create_atoms CPU = 0.000 seconds
group fluid type 1
300 atoms in group fluid
group rig type 2
100 atoms in group rig
displace_atoms fluid random 0.1 0.1 0 135414 units box
Displacing atoms ...
# ------ Model parameters ------ #
variable rho0 equal 1.0
variable cs equal 1.0
variable mp equal ${rho0}/${n}
variable mp equal 1/${n}
variable mp equal 1/1
variable zeta equal 1.0
variable kappa equal 1.0*${rho0}/${mp}
variable kappa equal 1.0*1/${mp}
variable kappa equal 1.0*1/1
variable fext equal 1e-4/${n}
variable fext equal 1e-4/1
variable dt_max equal 0.1*${cut}/${cs}/3
variable dt_max equal 0.1*3/${cs}/3
variable dt_max equal 0.1*3/1/3
variable Dr equal 0.05*${cut}*${cs}
variable Dr equal 0.05*3*${cs}
variable Dr equal 0.05*3*1
variable eta equal 0.1
variable gd0 equal 5e-4
variable npow equal 0.5
variable K equal 0.001
mass * ${mp}
mass * 1
set group all rheo/rho ${rho0}
set group all rheo/rho 1
Setting atom values ...
400 settings made for rheo/rho
set group all rheo/status 0
Setting atom values ...
400 settings made for rheo/status
set group rig rheo/status 1
Setting atom values ...
100 settings made for rheo/status
timestep ${dt_max}
timestep 0.1
pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
pair_style rheo 3 artificial/visc ${zeta} rho/damp ${Dr}
pair_style rheo 3 artificial/visc 1 rho/damp ${Dr}
pair_style rheo 3 artificial/visc 1 rho/damp 0.15
pair_coeff * *
# ------ Fixes & computes ------ #
fix 1 all rheo ${cut} quintic 0 shift
fix 1 all rheo 3 quintic 0 shift
fix 2 all rheo/viscosity * constant ${eta}
fix 2 all rheo/viscosity * constant 0.1
#fix 2 all rheo/viscosity * power ${eta} ${gd0} ${K} ${npow}
fix 3 all rheo/pressure * linear
fix 4 rig setforce 0.0 0.0 0.0
fix 5 fluid addforce ${fext} 0.0 0.0
fix 5 fluid addforce 0.0001 0.0 0.0
fix 6 all enforce2d
compute rho all rheo/property/atom rho
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press
dump 1 all custom 200 atomDump id type x y vx vy fx fy c_rho
run 20000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- @article{PalermoInPrep,
journal = {in prep},
title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},
year = {2024},
author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
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 = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 13 13 1
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair rheo, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) compute RHEO/KERNEL, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/2d
bin: standard
(3) compute RHEO/GRAD, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(4) compute RHEO/VSHIFT, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 5.693 | 5.693 | 5.693 Mbytes
Step Time KinEng Press
0 0 0 0
200 20 1.2220462e-06 3.7383146e-05
400 40 4.345762e-06 7.5866885e-05
600 60 8.8559433e-06 0.00011353743
800 80 1.4370506e-05 0.00015135634
1000 100 2.0576198e-05 0.00018903722
1200 120 2.721926e-05 0.00022533997
1400 140 3.4099653e-05 0.00026016069
1600 160 4.1064175e-05 0.00029445207
1800 180 4.8001225e-05 0.00032893763
2000 200 5.4832849e-05 0.00036402396
2200 220 6.1508431e-05 0.00039945249
2400 240 6.8000141e-05 0.00043534411
2600 260 7.430136e-05 0.00046943441
2800 280 8.0415328e-05 0.00049807225
3000 300 8.6335032e-05 0.00051815375
3200 320 9.2021626e-05 0.00052618224
3400 340 9.7387936e-05 0.00051877918
3600 360 0.00010231526 0.00048650828
3800 380 0.00010676617 0.00044578079
4000 400 0.00011080098 0.00044777126
4200 420 0.00011448127 0.00047047629
4400 440 0.00011787852 0.00050280249
4600 460 0.00012106805 0.0005397213
4800 480 0.00012412056 0.00057885539
5000 500 0.0001271078 0.00061396896
5200 520 0.00013006637 0.00063981812
5400 540 0.00013295039 0.00065094073
5600 560 0.00013561487 0.00063918847
5800 580 0.00013791796 0.00059087656
6000 600 0.00013983422 0.00052171998
6200 620 0.00014144833 0.00050658002
6400 640 0.00014286538 0.0005248626
6600 660 0.00014417734 0.00055826606
6800 680 0.00014546931 0.00060063748
7000 700 0.00014682553 0.00064421411
7200 720 0.0001482833 0.00068252242
7400 740 0.00014977996 0.00070671308
7600 760 0.00015114829 0.00069774026
7800 780 0.0001522719 0.00064408311
8000 800 0.00015312897 0.00055977044
8200 820 0.00015375669 0.0005225573
8400 840 0.00015425683 0.00053833691
8600 860 0.00015471278 0.00057447427
8800 880 0.0001552059 0.00061980921
9000 900 0.00015581593 0.0006659836
9200 920 0.0001565564 0.00070813532
9400 940 0.00015733573 0.00073378551
9600 960 0.00015802107 0.00071560835
9800 980 0.00015855339 0.00065636189
10000 1000 0.00015890743 0.0005699855
10200 1020 0.00015908095 0.00053138971
10400 1040 0.00015915523 0.00054790708
10600 1060 0.00015921254 0.00058899454
10800 1080 0.00015934193 0.00063964906
11000 1100 0.00015959891 0.00069241358
11200 1120 0.0001599636 0.00073734651
11400 1140 0.00016036526 0.00074477329
11600 1160 0.00016075471 0.00071047555
11800 1180 0.00016109516 0.00064173183
12000 1200 0.00016131524 0.00055500553
12200 1220 0.00016136366 0.0005290215
12400 1240 0.0001613025 0.00055124296
12600 1260 0.00016123023 0.00059758627
12800 1280 0.00016123043 0.00065488735
13000 1300 0.00016132935 0.0007140876
13200 1320 0.00016152165 0.00074795629
13400 1340 0.00016180372 0.00074730778
13600 1360 0.00016216585 0.00071370995
13800 1380 0.0001625339 0.00065176323
14000 1400 0.00016274999 0.00057515371
14200 1420 0.00016271295 0.00055878258
14400 1440 0.00016249768 0.00058448193
14600 1460 0.00016223675 0.00063096229
14800 1480 0.00016201846 0.00068639548
15000 1500 0.00016190593 0.00072444357
15200 1520 0.00016194466 0.00073830636
15400 1540 0.00016216164 0.00072773256
15600 1560 0.00016253174 0.00069215481
15800 1580 0.00016290895 0.00063239408
16000 1600 0.00016306463 0.00057466273
16200 1620 0.00016292218 0.00057951567
16400 1640 0.00016261117 0.00061504156
16600 1660 0.00016225906 0.00066066637
16800 1680 0.00016197993 0.00069751908
17000 1700 0.0001618568 0.00072202303
17200 1720 0.00016194264 0.00073255034
17400 1740 0.00016225911 0.0007231031
17600 1760 0.00016270465 0.00068931224
17800 1780 0.00016304053 0.00062934836
18000 1800 0.00016302624 0.00058060272
18200 1820 0.00016274847 0.00058859513
18400 1840 0.00016236893 0.00061804803
18600 1860 0.00016202777 0.00065393237
18800 1880 0.0001618184 0.00068747094
19000 1900 0.0001618044 0.00071352541
19200 1920 0.00016204402 0.00072351769
19400 1940 0.00016249999 0.00071330322
19600 1960 0.00016297924 0.00067984167
19800 1980 0.00016317435 0.00061634142
20000 2000 0.00016301186 0.00057234115
Loop time of 15.6198 on 4 procs for 20000 steps with 400 atoms
Performance: 11062881.511 tau/day, 1280.426 timesteps/s, 512.170 katom-step/s
99.7% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.1979 | 2.4473 | 2.6992 | 15.7 | 15.67
Neigh | 0.024709 | 0.027006 | 0.029223 | 1.3 | 0.17
Comm | 0.4657 | 0.71686 | 0.9662 | 29.0 | 4.59
Output | 0.033698 | 0.036781 | 0.039359 | 1.1 | 0.24
Modify | 12.306 | 12.313 | 12.319 | 0.2 | 78.83
Other | | 0.07916 | | | 0.51
Nlocal: 100 ave 107 max 93 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Nghost: 185.5 ave 192 max 179 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Neighs: 1712 ave 1848 max 1598 min
Histogram: 1 0 1 0 0 1 0 0 0 1
FullNghs: 3424 ave 3682 max 3174 min
Histogram: 1 0 1 0 0 0 1 0 0 1
Total # of neighbors = 13696
Ave neighs/atom = 34.24
Neighbor list builds = 331
Dangerous builds = 0
Total wall time: 0:00:15

65
examples/rheo/taylor-green/in.rheo.taylor.green Normal file
View File

@ -0,0 +1,65 @@
# ------ 2D Taylor Green vortex ------ #
dimension 2
units lj
atom_style rheo
boundary p p p
comm_modify vel yes
newton off
# ------ Create simulation box ------ #
variable n equal 1.0
variable cut equal 3.0
region box block 0 40 0 40 -0.01 0.01
create_box 1 box
lattice sq ${n}
create_atoms 1 region box
displace_atoms all random 0.1 0.1 0 135414 units box
# ------ Model parameters ------ #
variable rho0 equal 1.0
variable mp equal ${rho0}/${n}
variable cs equal 1.0
variable eta equal 0.05
variable zeta equal 1
variable dt_max equal 0.1*${cut}/${cs}/3
variable Dr equal 0.1*${cut}*${cs}
mass * ${mp}
set group all rheo/rho ${rho0}
set group all rheo/status 0
variable u0 equal 0.05
variable uy atom ${u0}*sin(2*PI*x/lx)*cos(2*PI*y/ly)
variable ux atom -${u0}*sin(2*PI*y/ly)*cos(2*PI*x/ly)
variable d0 atom ${rho0}-${u0}*${u0}*${rho0}*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/${cs}/${cs}
velocity all set v_ux v_uy 0.0 units box
timestep ${dt_max}
pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
pair_coeff * *
# ------ Fixes & computes ------ #
fix 1 all rheo ${cut} RK1 8 shift
fix 2 all rheo/viscosity * constant ${eta}
fix 3 all rheo/pressure * linear
fix 4 all enforce2d
compute rho all rheo/property/atom rho
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_rho
run 10000

224
examples/rheo/taylor-green/log.17Apr2024.taylor.green.g++.4 Normal file
View File

@ -0,0 +1,224 @@
LAMMPS (17 Apr 2024 - Development - patch_5May2020-18508-g3c0eaf6870-modified)
# ------ 2D Taylor Green vortex ------ #
dimension 2
units lj
atom_style rheo
boundary p p p
comm_modify vel yes
newton off
# ------ Create simulation box ------ #
variable n equal 1.0
variable cut equal 3.0
region box block 0 40 0 40 -0.01 0.01
create_box 1 box
Created orthogonal box = (0 0 -0.01) to (40 40 0.01)
2 by 2 by 1 MPI processor grid
lattice sq ${n}
lattice sq 1
Lattice spacing in x,y,z = 1 1 1
create_atoms 1 region box
Created 1600 atoms
using lattice units in orthogonal box = (0 0 -0.01) to (40 40 0.01)
create_atoms CPU = 0.001 seconds
displace_atoms all random 0.1 0.1 0 135414 units box
Displacing atoms ...
# ------ Model parameters ------ #
variable rho0 equal 1.0
variable mp equal ${rho0}/${n}
variable mp equal 1/${n}
variable mp equal 1/1
variable cs equal 1.0
variable eta equal 0.05
variable zeta equal 1
variable dt_max equal 0.1*${cut}/${cs}/3
variable dt_max equal 0.1*3/${cs}/3
variable dt_max equal 0.1*3/1/3
variable Dr equal 0.1*${cut}*${cs}
variable Dr equal 0.1*3*${cs}
variable Dr equal 0.1*3*1
mass * ${mp}
mass * 1
set group all rheo/rho ${rho0}
set group all rheo/rho 1
Setting atom values ...
1600 settings made for rheo/rho
set group all rheo/status 0
Setting atom values ...
1600 settings made for rheo/status
variable u0 equal 0.05
variable uy atom ${u0}*sin(2*PI*x/lx)*cos(2*PI*y/ly)
variable uy atom 0.05*sin(2*PI*x/lx)*cos(2*PI*y/ly)
variable ux atom -${u0}*sin(2*PI*y/ly)*cos(2*PI*x/ly)
variable ux atom -0.05*sin(2*PI*y/ly)*cos(2*PI*x/ly)
variable d0 atom ${rho0}-${u0}*${u0}*${rho0}*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/${cs}/${cs}
variable d0 atom 1-${u0}*${u0}*${rho0}*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/${cs}/${cs}
variable d0 atom 1-0.05*${u0}*${rho0}*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/${cs}/${cs}
variable d0 atom 1-0.05*0.05*${rho0}*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/${cs}/${cs}
variable d0 atom 1-0.05*0.05*1*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/${cs}/${cs}
variable d0 atom 1-0.05*0.05*1*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/1/${cs}
variable d0 atom 1-0.05*0.05*1*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/1/1
velocity all set v_ux v_uy 0.0 units box
timestep ${dt_max}
timestep 0.1
pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
pair_style rheo 3 artificial/visc ${zeta} rho/damp ${Dr}
pair_style rheo 3 artificial/visc 1 rho/damp ${Dr}
pair_style rheo 3 artificial/visc 1 rho/damp 0.3
pair_coeff * *
# ------ Fixes & computes ------ #
fix 1 all rheo ${cut} RK1 8 shift
fix 1 all rheo 3 RK1 8 shift
fix 2 all rheo/viscosity * constant ${eta}
fix 2 all rheo/viscosity * constant 0.05
fix 3 all rheo/pressure * linear
fix 4 all enforce2d
compute rho all rheo/property/atom rho
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press
dump 1 all custom 200 atomDump id type x y vx vy fx fy c_rho
run 10000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- @article{PalermoInPrep,
journal = {in prep},
title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},
year = {2024},
author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 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 = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 25 25 1
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair rheo, perpetual, half/full from (2)
attributes: half, newton off
pair build: halffull/newtoff
stencil: none
bin: none
(2) compute RHEO/KERNEL, perpetual
attributes: full, newton off
pair build: full/bin/atomonly
stencil: full/bin/2d
bin: standard
(3) compute RHEO/GRAD, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(4) compute RHEO/VSHIFT, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 6.835 | 6.835 | 6.835 Mbytes
Step Time KinEng Press
0 0 0.00062497276 0.00062607301
200 20 0.00056200647 0.00056633785
400 40 0.00050570968 0.00051098771
600 60 0.00045586684 0.00046081672
800 80 0.00041124523 0.00041549607
1000 100 0.00037065341 0.00037412741
1200 120 0.00033391585 0.00033580899
1400 140 0.00030078316 0.00030057307
1600 160 0.00027093231 0.00026842603
1800 180 0.00024403239 0.00023839026
2000 200 0.0002197865 0.00021148941
2200 220 0.0001979269 0.00018659386
2400 240 0.00017822267 0.00016430442
2600 260 0.00016047141 0.00014408514
2800 280 0.00014448504 0.00012574125
3000 300 0.00013009159 0.00010869938
3200 320 0.00011713578 9.414951e-05
3400 340 0.00010547564 8.1900579e-05
3600 360 9.4982139e-05 7.1285649e-05
3800 380 8.5538983e-05 6.1571123e-05
4000 400 7.7040171e-05 5.3462572e-05
4200 420 6.9390317e-05 4.6338308e-05
4400 440 6.2503763e-05 3.9697323e-05
4600 460 5.6303766e-05 3.4234465e-05
4800 480 5.0721595e-05 3.0841338e-05
5000 500 4.5695301e-05 2.7788566e-05
5200 520 4.1169161e-05 2.5744409e-05
5400 540 3.7093059e-05 2.3912739e-05
5600 560 3.3421819e-05 2.2494185e-05
5800 580 3.0114735e-05 2.1594384e-05
6000 600 2.7135224e-05 2.1164421e-05
6200 620 2.4450446e-05 2.0979349e-05
6400 640 2.2030925e-05 2.0858567e-05
6600 660 1.9850196e-05 2.098115e-05
6800 680 1.7884553e-05 2.1134827e-05
7000 700 1.6112763e-05 2.1242242e-05
7200 720 1.4515783e-05 2.1312763e-05
7400 740 1.3076456e-05 2.1370947e-05
7600 760 1.1779327e-05 2.1332126e-05
7800 780 1.0610469e-05 2.1156562e-05
8000 800 9.5573298e-06 2.0898126e-05
8200 820 8.6085799e-06 2.0517958e-05
8400 840 7.7539888e-06 1.9841551e-05
8600 860 6.9843033e-06 1.9114769e-05
8800 880 6.2911575e-06 1.8362959e-05
9000 900 5.6669785e-06 1.7473404e-05
9200 920 5.1049208e-06 1.6452745e-05
9400 940 4.5987908e-06 1.5578629e-05
9600 960 4.1429972e-06 1.4427274e-05
9800 980 3.7324962e-06 1.3169649e-05
10000 1000 3.3627455e-06 1.1938723e-05
Loop time of 38.2006 on 4 procs for 10000 steps with 1600 atoms
Performance: 2261743.875 tau/day, 261.776 timesteps/s, 418.841 katom-step/s
99.7% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 8.2958 | 8.7273 | 9.3582 | 15.2 | 22.85
Neigh | 0.034282 | 0.035689 | 0.037115 | 0.7 | 0.09
Comm | 0.16788 | 0.17018 | 0.17278 | 0.4 | 0.45
Output | 0.066977 | 0.06882 | 0.071704 | 0.7 | 0.18
Modify | 28.483 | 28.793 | 28.962 | 3.6 | 75.37
Other | | 0.4053 | | | 1.06
Nlocal: 400 ave 402 max 399 min
Histogram: 2 0 0 1 0 0 0 0 0 1
Nghost: 307.25 ave 308 max 305 min
Histogram: 1 0 0 0 0 0 0 0 0 3
Neighs: 7618.25 ave 7697 max 7564 min
Histogram: 1 0 1 1 0 0 0 0 0 1
FullNghs: 13343 ave 13497 max 13258 min
Histogram: 1 1 1 0 0 0 0 0 0 1
Total # of neighbors = 53372
Ave neighs/atom = 33.3575
Neighbor list builds = 123
Dangerous builds = 0
Total wall time: 0:00:38

4
lib/qmmm/Makefile.gfortran-cmake
View File

@ -46,8 +46,8 @@ QELIBS += -lgfortran -lmpi_mpifh
# part 3: add-on libraries and main library for LAMMPS
sinclude ../../src/Makefile.package
LAMMPSFLAGS = $(shell pkgconf liblammps --cflags)
LAMMPSLIB = $(shell pkgconf liblammps --libs)
LAMMPSFLAGS = $(shell pkg-config --cflags liblammps)
LAMMPSLIB = $(shell pkg-config --libs liblammps)
# part 4: local QM/MM library and progams
SRC=pwqmmm.c libqmmm.c

2
lib/qmmm/README
View File

@ -101,7 +101,7 @@ Makefile.gfortran-cmake and make adjustments to the makefile variables
according to the comments in the file. You probably need to adjust
the QETOPDIR variable to point to the location of your QE
compilation/installation.
Please also check that the command "pkgconf liblammps --libs" works.
Please also check that the command "pkg-config --libs liblammps" works.
Then you should be able to compile the QM/MM executable with:
make -f Makefile.gfortran-cmake pwqmmm.x

14
lib/rheo/Makefile.lammps Normal file
View File

@ -0,0 +1,14 @@
# Settings that the LAMMPS build will import when this package is installed
ifeq ($(strip $(shell pkg-config --version)),)
# manual configuration w/o pkg-config/pkgconf
# change this to -I/path/to/your/lib/gsl/include/
rheo_SYSINC = -I../../lib/rheo/gsl/include/
# change this to -L/path/to/your/lib/gsl/lib/
rheo_SYSLIB = -L../../lib/rheo/gsl/lib/ -lgsl -lgslcblas
else
# autodetect GSL settings from pkg-config/pkgconf
rheo_SYSINC = $(shell pkg-config --cflags gsl)
rheo_SYSLIB = $(shell pkg-config --libs gsl)
endif

7
lib/rheo/README Normal file
View File

@ -0,0 +1,7 @@
This directory has a Makefile.lammps file with settings that allows LAMMPS to
dynamically link to the GSL library. This is required to use the RHEO package
in a LAMMPS input script. If you have the pkg-config command available, it
will automatically import the GSL settings. Otherwise they will have to be
added manually.
See the header of Makefile.lammps for more info.

5
src/.gitignore vendored
View File

@ -245,6 +245,9 @@
/pair_tdpd.cpp
/pair_tdpd.h
/*rheo*.cpp
/*rheo*.h
/compute_grid.cpp
/compute_grid.h
/compute_grid_local.cpp
@ -790,6 +793,8 @@
/fix_acks2_reaxff.h
/fix_adapt_fep.cpp
/fix_adapt_fep.h
/fix_add_heat.cpp
/fix_add_heat.h
/fix_addtorque.cpp
/fix_addtorque.h
/fix_append_atoms.cpp

58
src/BPM/bond_bpm.cpp
View File

@ -18,6 +18,7 @@
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "fix_bond_history.h"
#include "fix_store_local.h"
#include "fix_update_special_bonds.h"
@ -53,10 +54,14 @@ BondBPM::BondBPM(LAMMPS *_lmp) :
pack_choice(nullptr), output_data(nullptr)
{
overlay_flag = 0;
ignore_special_flag = 0;
prop_atom_flag = 0;
break_flag = 1;
nvalues = 0;
nhistory = 0;
update_flag = 0;
r0_max_estimate = 0.0;
max_stretch = 1.0;
@ -64,10 +69,10 @@ BondBPM::BondBPM(LAMMPS *_lmp) :
// this is so final order of Modify:fix will conform to input script
// BondHistory technically only needs this if updateflag = 1
id_fix_dummy = utils::strdup("BPM_DUMMY");
id_fix_dummy = utils::strdup(fmt::format("BPM_DUMMY_{}", instance_total));
modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy));
id_fix_dummy2 = utils::strdup("BPM_DUMMY2");
id_fix_dummy2 = utils::strdup(fmt::format("BPM_DUMMY2_{}", instance_total));
modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy2));
if (lmp->citeme) lmp->citeme->add(cite_bpm);
@ -82,7 +87,7 @@ BondBPM::~BondBPM()
if (id_fix_dummy) modify->delete_fix(id_fix_dummy);
if (id_fix_dummy2) modify->delete_fix(id_fix_dummy2);
if (id_fix_update) modify->delete_fix(id_fix_update);
if (id_fix_bond_history) modify->delete_fix(id_fix_bond_history);
if (fix_bond_history) modify->delete_fix(id_fix_bond_history);
if (id_fix_store_local) modify->delete_fix(id_fix_store_local);
if (id_fix_prop_atom) modify->delete_fix(id_fix_prop_atom);
@ -109,6 +114,7 @@ void BondBPM::init_style()
fix_store_local->nvalues = nvalues;
}
if (!ignore_special_flag) {
if (overlay_flag) {
if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 ||
force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0)
@ -144,15 +150,27 @@ void BondBPM::init_style()
}
}
// special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists
if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || force->special_coul[2] != 1.0 ||
force->special_coul[3] != 1.0)
error->all(FLERR, "Bond style bpm requires 1-3 and 1-4 special weights of 1.0");
}
if (force->angle || force->dihedral || force->improper)
error->all(FLERR, "Bond style bpm cannot be used with 3,4-body interactions");
if (atom->molecular == 2)
error->all(FLERR, "Bond style bpm cannot be used with atom style template");
// special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists
if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || force->special_coul[2] != 1.0 ||
force->special_coul[3] != 1.0)
error->all(FLERR, "Bond style bpm requires 1-3 and 1-4 special weights of 1.0");
// find all instances of bond history to delete/shift data
// (bond hybrid may create multiple)
histories = modify->get_fix_by_style("BOND_HISTORY");
n_histories = histories.size();
// If a bond type isn't set, must be using bond style hybrid
hybrid_flag = 0;
for (int i = 1; i <= atom->nbondtypes; i++)
if (!setflag[i]) hybrid_flag = 1;
fix_bond_history->setflag = setflag;
}
/* ----------------------------------------------------------------------
@ -269,6 +287,14 @@ void BondBPM::settings(int narg, char **arg)
}
}
}
// Set up necessary history fix
if (!fix_bond_history) {
fix_bond_history = dynamic_cast<FixBondHistory *>(modify->replace_fix(
id_fix_dummy2, fmt::format("{} all BOND_HISTORY {} {}", id_fix_bond_history, update_flag, nhistory), 1));
delete[] id_fix_dummy2;
id_fix_dummy2 = nullptr;
}
}
/* ----------------------------------------------------------------------
@ -385,12 +411,15 @@ void BondBPM::process_broken(int i, int j)
if (i < nlocal) {
for (m = 0; m < num_bond[i]; m++) {
if (bond_atom[i][m] == tag[j]) {
if (bond_atom[i][m] == tag[j] && setflag[bond_type[i][m]]) {
n = num_bond[i];
bond_type[i][m] = bond_type[i][n - 1];
bond_atom[i][m] = bond_atom[i][n - 1];
fix_bond_history->shift_history(i, m, n - 1);
fix_bond_history->delete_history(i, n - 1);
for (auto &ihistory: histories) {
auto fix_bond_history2 = dynamic_cast<FixBondHistory *> (ihistory);
fix_bond_history2->shift_history(i, m, n - 1);
fix_bond_history2->delete_history(i, n - 1);
}
num_bond[i]--;
break;
}
@ -399,12 +428,15 @@ void BondBPM::process_broken(int i, int j)
if (j < nlocal) {
for (m = 0; m < num_bond[j]; m++) {
if (bond_atom[j][m] == tag[i]) {
if (bond_atom[j][m] == tag[i] && setflag[bond_type[j][m]]) {
n = num_bond[j];
bond_type[j][m] = bond_type[j][n - 1];
bond_atom[j][m] = bond_atom[j][n - 1];
fix_bond_history->shift_history(j, m, n - 1);
fix_bond_history->delete_history(j, n - 1);
for (auto &ihistory: histories) {
auto fix_bond_history2 = dynamic_cast<FixBondHistory *> (ihistory);
fix_bond_history2->shift_history(j, m, n - 1);
fix_bond_history2->delete_history(j, n - 1);
}
num_bond[j]--;
break;
}

11
src/BPM/bond_bpm.h
View File

@ -16,8 +16,12 @@
#include "bond.h"
#include <vector>
namespace LAMMPS_NS {
class Fix;
class BondBPM : public Bond {
public:
BondBPM(class LAMMPS *);
@ -34,7 +38,7 @@ class BondBPM : public Bond {
protected:
double r0_max_estimate;
double max_stretch;
int store_local_freq;
int store_local_freq, nhistory, update_flag, hybrid_flag;
std::vector<int> leftover_iarg;
@ -50,9 +54,12 @@ class BondBPM : public Bond {
FnPtrPack *pack_choice; // ptrs to pack functions
double *output_data;
int prop_atom_flag, nvalues, overlay_flag, break_flag;
int prop_atom_flag, nvalues, overlay_flag, break_flag, ignore_special_flag;
int index_x_ref, index_y_ref, index_z_ref;
int n_histories;
std::vector<Fix *> histories;
void pack_id1(int, int, int);
void pack_id2(int, int, int);
void pack_time(int, int, int);

17
src/BPM/bond_bpm_rotational.cpp
View File

@ -52,6 +52,9 @@ BondBPMRotational::BondBPMRotational(LAMMPS *_lmp) :
smooth_flag = 1;
normalize_flag = 0;
nhistory = 4;
id_fix_bond_history = utils::strdup("HISTORY_BPM_ROTATIONAL");
single_extra = 7;
svector = new double[7];
}
@ -458,6 +461,9 @@ void BondBPMRotational::compute(int eflag, int vflag)
store_data();
}
if (hybrid_flag)
fix_bond_history->compress_history();
int i1, i2, itmp, n, type;
double r[3], r0[3], rhat[3];
double rsq, r0_mag, r_mag, r_mag_inv;
@ -563,6 +569,9 @@ void BondBPMRotational::compute(int eflag, int vflag)
ev_tally_xyz(i1, i2, nlocal, newton_bond, 0.0, -force1on2[0] * smooth, -force1on2[1] * smooth,
-force1on2[2] * smooth, r[0], r[1], r[2]);
}
if (hybrid_flag)
fix_bond_history->uncompress_history();
}
/* ---------------------------------------------------------------------- */
@ -652,14 +661,6 @@ void BondBPMRotational::init_style()
if (domain->dimension == 2)
error->warning(FLERR, "Bond style bpm/rotational not intended for 2d use");
if (!id_fix_bond_history) {
id_fix_bond_history = utils::strdup("HISTORY_BPM_ROTATIONAL");
fix_bond_history = dynamic_cast<FixBondHistory *>(modify->replace_fix(
id_fix_dummy2, fmt::format("{} all BOND_HISTORY 0 4", id_fix_bond_history), 1));
delete[] id_fix_dummy2;
id_fix_dummy2 = nullptr;
}
}
/* ---------------------------------------------------------------------- */

17
src/BPM/bond_bpm_spring.cpp
View File

@ -39,6 +39,9 @@ BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) :
smooth_flag = 1;
normalize_flag = 0;
nhistory = 1;
id_fix_bond_history = utils::strdup("HISTORY_BPM_SPRING");
single_extra = 1;
svector = new double[1];
}
@ -137,6 +140,9 @@ void BondBPMSpring::compute(int eflag, int vflag)
store_data();
}
if (hybrid_flag)
fix_bond_history->compress_history();
int i1, i2, itmp, n, type;
double delx, dely, delz, delvx, delvy, delvz;
double e, rsq, r, r0, rinv, smooth, fbond, dot;
@ -226,6 +232,9 @@ void BondBPMSpring::compute(int eflag, int vflag)
if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, fbond, delx, dely, delz);
}
if (hybrid_flag)
fix_bond_history->uncompress_history();
}
/* ---------------------------------------------------------------------- */
@ -283,14 +292,6 @@ void BondBPMSpring::init_style()
if (comm->ghost_velocity == 0)
error->all(FLERR, "Bond bpm/spring requires ghost atoms store velocity");
if (!id_fix_bond_history) {
id_fix_bond_history = utils::strdup("HISTORY_BPM_SPRING");
fix_bond_history = dynamic_cast<FixBondHistory *>(modify->replace_fix(
id_fix_dummy2, fmt::format("{} all BOND_HISTORY 0 1", id_fix_bond_history), 1));
delete[] id_fix_dummy2;
id_fix_dummy2 = nullptr;
}
}
/* ---------------------------------------------------------------------- */

18
src/BPM/compute_nbond_atom.cpp
View File

@ -14,7 +14,9 @@
#include "compute_nbond_atom.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
@ -27,6 +29,21 @@ ComputeNBondAtom::ComputeNBondAtom(LAMMPS *_lmp, int narg, char **arg) :
{
if (narg < 3) utils::missing_cmd_args(FLERR, "compute nbond/atom", error);
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Compute nbond/atom used in system without bonds");
btype = -1;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg], "bond/type") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute nbond/atom bond/type", error);
btype = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else {
error->all(FLERR, "Unknown compute nbond/atom command {}", arg[iarg]);
}
}
peratom_flag = 1;
size_peratom_cols = 0;
comm_reverse = 1;
@ -77,6 +94,7 @@ void ComputeNBondAtom::compute_peratom()
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_bond[i]; j++) {
if (bond_type[i][j] <= 0) continue;
if (btype != -1 && bond_type[i][j] != btype) continue;
k = atom->map(bond_atom[i][j]);
if (k < 0) continue;

2
src/BPM/compute_nbond_atom.h
View File

@ -35,7 +35,7 @@ class ComputeNBondAtom : public Compute {
double memory_usage() override;
private:
int nmax;
int nmax, btype;
double *nbond;
};

13
src/BPM/fix_update_special_bonds.cpp
View File

@ -90,11 +90,13 @@ void FixUpdateSpecialBonds::pre_exchange()
for (auto const &it : broken_pairs) {
tagi = it.first;
tagj = it.second;
i = atom->map(tagi);
j = atom->map(tagj);
// remove i from special bond list for atom j and vice versa
// ignore n2, n3 since 1-3, 1-4 special factors required to be 1.0
// assume ghosts don't need special information
if (i < nlocal) {
slist = special[i];
n1 = nspecial[i][0];
@ -126,20 +128,25 @@ void FixUpdateSpecialBonds::pre_exchange()
// add i to special bond list for atom j and vice versa
// ignore n2, n3 since 1-3, 1-4 special factors required to be 1.0
// assume ghosts don't need special information
if (i < nlocal) {
n1 = nspecial[i][0];
if (n1 >= atom->maxspecial)
error->one(FLERR, "Special list size exceeded in fix update/special/bond");
error->one(FLERR, "Special list size exceeded for atom {}", tagi);
special[i][n1] = tagj;
nspecial[i][0] += 1;
nspecial[i][1] = nspecial[i][2] = nspecial[i][0];
}
if (j < nlocal) {
n1 = nspecial[j][0];
if (n1 >= atom->maxspecial)
error->one(FLERR, "Special list size exceeded in fix update/special/bond");
error->one(FLERR, "Special list size exceeded for atom {}", tagj);
special[j][n1] = tagi;
nspecial[j][0] += 1;
nspecial[j][1] = nspecial[j][2] = nspecial[j][0];
}
}
broken_pairs.clear();
created_pairs.clear();
@ -162,7 +169,7 @@ void FixUpdateSpecialBonds::pre_force(int /*vflag*/)
tagint *tag = atom->tag;
// In theory could communicate a list of broken bonds to neighboring processors here
// to remove restriction that users use Newton bond off
// to remove restriction on Newton bond off
for (int ilist = 0; ilist < neighbor->nlist; ilist++) {
list = neighbor->lists[ilist];

131
src/ELECTRODE/fix_electrode_conp.cpp
View File

@ -79,6 +79,9 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
potential_i(nullptr), potential_iele(nullptr)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_electrode);
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR, "Fix {} requires an atom map, see atom_modify", style);
// fix.h output flags
scalar_flag = 1;
vector_flag = 1;
@ -86,6 +89,9 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
extvector = 0;
extarray = 0;
virial_global_flag = 1; // use virials of this fix
thermo_virial = 1; // set vflags for v_tally
bool default_algo = true;
algo = Algo::MATRIX_INV;
matrix_algo = true;
@ -123,7 +129,8 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
} else
group_psi_const[0] = utils::numeric(FLERR, arg[3], false, lmp);
char *eta_str = arg[4];
eta = utils::numeric(FLERR, eta_str, false, lmp);
bool etanull = (strcmp(eta_str, "NULL") == 0);
if (!etanull) eta = utils::numeric(FLERR, eta_str, false, lmp);
int iarg = 5;
while (iarg < narg) {
if ((strcmp(arg[iarg], "couple") == 0)) {
@ -214,7 +221,7 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
qtotal_var_style = VarStyle::CONST;
}
} else if ((strcmp(arg[iarg], "eta") == 0)) {
if (iarg + 2 > narg) error->all(FLERR, "Need two arguments after eta command");
if (iarg + 2 > narg) error->all(FLERR, "Need one argument after eta command");
etaflag = true;
int is_double, cols, ghost;
eta_index = atom->find_custom_ghost(arg[++iarg] + 2, is_double, cols, ghost);
@ -240,6 +247,7 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
if (qtotal_var_style != VarStyle::UNSET) {
if (symm) error->all(FLERR, "{} cannot use qtotal keyword with symm on", this->style);
}
if (etanull && !etaflag) error->all(FLERR, "If eta is NULL the eta keyword must be used");
// computatonal potential
group_psi = std::vector<double>(groups.size());
@ -346,7 +354,7 @@ int FixElectrodeConp::modify_param(int narg, char **arg)
MPI_Allreduce(MPI_IN_PLACE, &in_ele, 1, MPI_INT, MPI_SUM, world);
if (in_ele == 0) error->all(FLERR, "No atoms of type in electrode");
MPI_Allreduce(MPI_IN_PLACE, &not_in_ele, 1, MPI_INT, MPI_SUM, world);
if (not_in_ele)
if (not_in_ele && (comm->me == 0))
error->warning(FLERR,
"Not all atoms of type in electrode; Thomas-Fermi parameters will be ignored "
"for electrolyte");
@ -427,11 +435,10 @@ void FixElectrodeConp::init()
}
if (comm->me == 0)
for (char *fix_id : integrate_ids)
error->warning(
FLERR,
"Electrode atoms are integrated by fix {}, but fix electrode is using a matrix method. "
"For "
"mobile electrodes use the conjugate gradient algorithm without matrix ('algo cg').",
error->warning(FLERR,
"Electrode atoms are integrated by fix {}, but fix electrode is using a "
"matrix method. For mobile electrodes use the conjugate gradient algorithm "
"without matrix ('algo cg').",
fix_id);
}
@ -478,7 +485,7 @@ void FixElectrodeConp::post_constructor()
input->variable->set(fmt::format("{} equal f_{}[{}]", var_vtop, fixname, 1 + top_group));
input->variable->set(fmt::format("{} equal (v_{}-v_{})/lz", var_efield, var_vbot, var_vtop));
// check for other efields and warn if found
if (modify->get_fix_by_style("^efield").size() > 0 && comm->me == 0)
if ((modify->get_fix_by_style("^efield").size() > 0) && (comm->me == 0))
error->warning(FLERR, "Other efield fixes found -- please make sure this is intended!");
// call fix command:
// fix [varstem]_efield all efield 0.0 0.0 [var_vdiff]/lz
@ -638,13 +645,14 @@ void FixElectrodeConp::setup_post_neighbor()
/* ---------------------------------------------------------------------- */
void FixElectrodeConp::setup_pre_reverse(int eflag, int /*vflag*/)
void FixElectrodeConp::setup_pre_reverse(int eflag, int vflag)
{
if (pair->did_tally_callback() && (comm->me == 0))
error->warning(FLERR, "Computation of virials in fix {} is incompatible with TALLY package", style);
// correct forces for initial timestep
gausscorr(eflag, true);
ev_init(eflag, vflag);
gausscorr(eflag, vflag, true);
self_energy(eflag);
// potential_energy(eflag); // not always part of the energy, depending on ensemble, therefore
// removed
}
/* ---------------------------------------------------------------------- */
@ -674,7 +682,7 @@ void FixElectrodeConp::invert()
void FixElectrodeConp::symmetrize()
{
// S matrix to enforce charge neutrality constraint
if (read_inv && comm->me == 0)
if (read_inv && (comm->me == 0))
error->warning(FLERR,
"Symmetrizing matrix from file. Make sure the provided matrix has not been "
"symmetrized yet.");
@ -755,11 +763,11 @@ void FixElectrodeConp::setup_pre_exchange() // create_taglist
// if memory_usage > 0.5 GiB, warn with expected usage
double mem_needed = memory_usage();
mem_needed /= (1024 * 1024 * 1024); // convert to GiB
if (mem_needed > 0.5 && comm->me == 0)
if ((mem_needed > 0.5) && (comm->me == 0))
error->warning(FLERR,
"Please ensure there is sufficient memory for fix electrode "
"Please ensure there is sufficient memory for fix {} "
"(anticipated usage is at least {:.1f} GiB per proc)",
mem_needed);
style, mem_needed);
}
/* ---------------------------------------------------------------------- */
@ -771,12 +779,11 @@ void FixElectrodeConp::pre_force(int)
/* ---------------------------------------------------------------------- */
void FixElectrodeConp::pre_reverse(int eflag, int /*vflag*/)
void FixElectrodeConp::pre_reverse(int eflag, int vflag)
{
gausscorr(eflag, true);
ev_init(eflag, vflag);
gausscorr(eflag, vflag, true);
self_energy(eflag);
//potential_energy(eflag); // not always part of the energy, depending on ensemble, therefore
// removed
}
/* ---------------------------------------------------------------------- */
@ -924,7 +931,7 @@ void FixElectrodeConp::update_charges()
dot_old = dot_new;
}
recompute_potential(std::move(b), q_local);
if (delta > cg_threshold && comm->me == 0) error->warning(FLERR, "CG threshold not reached");
if ((delta > cg_threshold) && (comm->me == 0)) error->warning(FLERR, "CG threshold not reached");
} else {
error->all(FLERR, "This algorithm is not implemented, yet");
}
@ -1224,11 +1231,10 @@ double FixElectrodeConp::self_energy(int eflag)
/* ---------------------------------------------------------------------- */
double FixElectrodeConp::gausscorr(int eflag, bool fflag)
double FixElectrodeConp::gausscorr(int eflag, int vflag, bool fflag)
{
// correction to short range interaction due to eta
int evflag = pair->evflag;
double const qqrd2e = force->qqrd2e;
int const nlocal = atom->nlocal;
int *mask = atom->mask;
@ -1294,13 +1300,11 @@ double FixElectrodeConp::gausscorr(int eflag, bool fflag)
f[j][2] -= delz * fpair;
}
}
double ecoul = 0.;
if (eflag) ecoul = -prefactor * erfc_etar;
if (evflag) {
force->pair->ev_tally(i, j, nlocal, newton_pair, 0., ecoul, fpair, delx, dely, delz);
if (eflag) {
double ecoul = -prefactor * erfc_etar;
force->pair->ev_tally(i, j, nlocal, newton_pair, 0., ecoul, 0., 0., 0., 0.);
}
if (vflag) v_tally(i, j, nlocal, newton_pair, fpair, delx, dely, delz);
}
}
}
@ -1621,3 +1625,70 @@ void FixElectrodeConp::unpack_forward_comm(int n, int first, double *buf)
int const last = first + n;
for (int i = first, m = 0; i < last; i++) atom->q[i] = buf[m++];
}
/* ----------------------------------------------------------------------
Tally virial of pair interactions in pre_reverse. This cannot be done with pair->ev_tally()
because compute_fdotr is called before pre_reverse, i.e. Virials need to be tallied even if fdotr
is used.
------------------------------------------------------------------------- */
void FixElectrodeConp::v_tally(int i, int j, int nlocal, int newton_pair, double fpair, double delx,
double dely, double delz)
{
double v[6];
if (vflag_either) {
v[0] = delx * delx * fpair;
v[1] = dely * dely * fpair;
v[2] = delz * delz * fpair;
v[3] = delx * dely * fpair;
v[4] = delx * delz * fpair;
v[5] = dely * delz * fpair;
if (vflag_global) {
if (newton_pair) {
virial[0] += v[0];
virial[1] += v[1];
virial[2] += v[2];
virial[3] += v[3];
virial[4] += v[4];
virial[5] += v[5];
} else {
if (i < nlocal) {
virial[0] += 0.5 * v[0];
virial[1] += 0.5 * v[1];
virial[2] += 0.5 * v[2];
virial[3] += 0.5 * v[3];
virial[4] += 0.5 * v[4];
virial[5] += 0.5 * v[5];
}
if (j < nlocal) {
virial[0] += 0.5 * v[0];
virial[1] += 0.5 * v[1];
virial[2] += 0.5 * v[2];
virial[3] += 0.5 * v[3];
virial[4] += 0.5 * v[4];
virial[5] += 0.5 * v[5];
}
}
}
if (vflag_atom) {
if (newton_pair || i < nlocal) {
vatom[i][0] += 0.5 * v[0];
vatom[i][1] += 0.5 * v[1];
vatom[i][2] += 0.5 * v[2];
vatom[i][3] += 0.5 * v[3];
vatom[i][4] += 0.5 * v[4];
vatom[i][5] += 0.5 * v[5];
}
if (newton_pair || j < nlocal) {
vatom[j][0] += 0.5 * v[0];
vatom[j][1] += 0.5 * v[1];
vatom[j][2] += 0.5 * v[2];
vatom[j][3] += 0.5 * v[3];
vatom[j][4] += 0.5 * v[4];
vatom[j][5] += 0.5 * v[5];
}
}
}
}

3
src/ELECTRODE/fix_electrode_conp.h
View File

@ -119,10 +119,11 @@ class FixElectrodeConp : public Fix {
void create_taglist();
void invert();
void symmetrize();
double gausscorr(int, bool);
double gausscorr(int, int, bool);
void update_charges();
double potential_energy();
double self_energy(int);
void v_tally(int, int, int, int, double, double, double, double);
void write_to_file(FILE *, const std::vector<tagint> &, const std::vector<std::vector<double>> &);
void read_from_file(const std::string &input_file, double **, const std::string &);
void compute_sd_vectors();

161
src/GRANULAR/fix_add_heat.cpp Normal file
View File

@ -0,0 +1,161 @@
/* ----------------------------------------------------------------------
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: Joel Clemmer (SNL)
------------------------------------------------------------------------- */
#include "fix_add_heat.h"
#include "atom.h"
#include "error.h"
#include "input.h"
#include "memory.h"
#include "update.h"
#include "variable.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum { CONSTANT, EQUAL, ATOM };
enum { ADD, LINEAR, QUARTIC };
/* ---------------------------------------------------------------------- */
FixAddHeat::FixAddHeat(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), varstr(nullptr), vatom(nullptr)
{
if (narg < 5) utils::missing_cmd_args(FLERR, "fix add/heat", error);
dynamic_group_allow = 1;
overwrite_flag = 0;
if (strcmp(arg[3], "constant") == 0) {
style = ADD;
} else if (strcmp(arg[3], "linear") == 0) {
style = LINEAR;
} else if (strcmp(arg[3], "quartic") == 0) {
style = QUARTIC;
} else {
error->all(FLERR, "Invalid option {}", arg[3]);
}
if (utils::strmatch(arg[4], "^v_")) {
varstr = utils::strdup(arg[4] + 2);
} else {
value = utils::numeric(FLERR, arg[4], false, lmp);
vstyle = CONSTANT;
}
int iarg = 5;
if (style != ADD) {
if (narg != 6) utils::missing_cmd_args(FLERR, "fix add/heat", error);
prefactor = utils::numeric(FLERR, arg[5], false, lmp);
iarg = 6;
}
// optional args
while (iarg < narg) {
if (strcmp(arg[iarg], "overwrite") == 0) {
if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix add/heat", error);
overwrite_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else {
error->all(FLERR, "Illegal fix add/heat command, invalid argument {}", arg[iarg]);
}
}
maxatom = -1;
}
/* ---------------------------------------------------------------------- */
FixAddHeat::~FixAddHeat()
{
delete[] varstr;
memory->destroy(vatom);
}
/* ---------------------------------------------------------------------- */
int FixAddHeat::setmask()
{
int mask = 0;
mask |= POST_FORCE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixAddHeat::init()
{
if (!atom->temperature_flag)
error->all(FLERR, "Fix add/heat requires atom style with temperature property");
if (!atom->heatflow_flag)
error->all(FLERR, "Fix add/heat requires atom style with heatflow property");
// check variable
if (varstr) {
var = input->variable->find(varstr);
if (var < 0) error->all(FLERR, "Variable {} for fix addforce does not exist", varstr);
if (input->variable->equalstyle(var))
vstyle = EQUAL;
else if (input->variable->atomstyle(var))
vstyle = ATOM;
else
error->all(FLERR, "Variable {} for fix addforce is invalid style", varstr);
}
}
/* ---------------------------------------------------------------------- */
void FixAddHeat::post_force(int /*vflag*/)
{
int *mask = atom->mask;
double *heatflow = atom->heatflow;
double *temperature = atom->temperature;
double dtinv = 1.0 / update->dt;
if (vstyle == ATOM) {
if (atom->nmax > maxatom) {
maxatom = atom->nmax;
memory->destroy(vatom);
memory->create(vatom, maxatom, "addheat:vatom");
}
input->variable->compute_atom(var, igroup, &vatom[0], 1, 0);
}
if (overwrite_flag)
for (int i = 0; i < atom->nlocal; i++)
if (mask[i] & groupbit)
heatflow[i] = 0.0;
double vtmp, dt;
if (vstyle == CONSTANT) vtmp = value;
if (vstyle == EQUAL) vtmp = input->variable->compute_equal(var);
for (int i = 0; i < atom->nlocal; i++) {
if (mask[i] & groupbit) {
if (vstyle == ATOM) vtmp = vatom[i];
if (style == ADD) {
heatflow[i] += dtinv * vtmp;
} else if (style == LINEAR) {
heatflow[i] += dtinv * prefactor * (vtmp - temperature[i]);
} else if (style == QUARTIC) {
heatflow[i] += dtinv * prefactor * (pow(vtmp, 4.0) - pow(temperature[i], 4.0));
}
}
}
}

45
src/GRANULAR/fix_add_heat.h Normal file
View File

@ -0,0 +1,45 @@
/* -*- 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 FIX_CLASS
// clang-format off
FixStyle(add/heat,FixAddHeat);
// clang-format on
#else
#ifndef LMP_FIX_ADD_HEAT_H
#define LMP_FIX_ADD_HEAT_H
#include "fix.h"
namespace LAMMPS_NS {
class FixAddHeat : public Fix {
public:
FixAddHeat(class LAMMPS *, int, char **);
~FixAddHeat() override;
int setmask() override;
void init() override;
void post_force(int) override;
protected:
double value, prefactor;
int var, vstyle, maxatom, style, overwrite_flag;
char *varstr;
double *vatom;
};
} // namespace LAMMPS_NS
#endif
#endif

3
src/Makefile
View File

@ -124,6 +124,7 @@ PACKAGE = \
reaction \
reaxff \
replica \
rheo \
rigid \
scafacos \
shock \
@ -233,6 +234,7 @@ PACKLIB = \
plumed \
qmmm \
ml-quip \
rheo \
scafacos \
machdyn \
vtk \
@ -254,6 +256,7 @@ PACKEXT = \
netcdf \
plumed \
qmmm \
rheo \
scafacos \
voronoi \
vtk \

7
src/PYTHON/fix_python_invoke.cpp
View File

@ -105,6 +105,13 @@ void FixPythonInvoke::end_of_step()
/* ---------------------------------------------------------------------- */
void FixPythonInvoke::setup(int vflag)
{
if (selected_callback == POST_FORCE) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixPythonInvoke::post_force(int vflag)
{
if (update->ntimestep % nevery != 0) return;

1
src/PYTHON/fix_python_invoke.h
View File

@ -30,6 +30,7 @@ class FixPythonInvoke : public Fix {
FixPythonInvoke(class LAMMPS *, int, char **);
~FixPythonInvoke() override;
int setmask() override;
void setup(int) override;
void end_of_step() override;
void post_force(int) override;

6
src/REAXFF/fix_reaxff_bonds.cpp
View File

@ -44,6 +44,7 @@ FixReaxFFBonds::FixReaxFFBonds(LAMMPS *lmp, int narg, char **arg) :
ntypes = atom->ntypes;
nmax = atom->nmax;
compressed = 0;
first_flag = true;
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
@ -94,7 +95,10 @@ int FixReaxFFBonds::setmask()
void FixReaxFFBonds::setup(int /*vflag*/)
{
end_of_step();
// only print output during setup() at the very beginning
// to avoid duplicate outputs when using multiple run statements
if (first_flag) end_of_step();
first_flag = false;
}
/* ---------------------------------------------------------------------- */

2
src/REAXFF/fix_reaxff_bonds.h
View File

@ -40,6 +40,7 @@ class FixReaxFFBonds : public Fix {
tagint **neighid;
double **abo;
FILE *fp;
bool first_flag;
void allocate();
void destroy();
@ -50,7 +51,6 @@ class FixReaxFFBonds : public Fix {
int nint(const double &);
double memory_usage() override;
bigint nvalid, nextvalid();
struct _reax_list *lists;
class PairReaxFF *reaxff;
class NeighList *list;

70
src/RHEO/Install.sh Normal file
View File

@ -0,0 +1,70 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# enforce using portable C locale
LC_ALL=C
export LC_ALL
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# some styles in RHEO have base classes in BPM
if (test $1 = 1) then
if (test ! -e ../bond_bpm.cpp) then
echo "Must install BPM package with RHEO"
exit 1
fi
fi
for file in *.cpp *.h; do
action ${file}
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*rheo[^ \t]* //' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(rheo_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(rheo_SYSLIB) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^[ \t]*include.*rheo.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/rheo\/Makefile.lammps
' ../Makefile.package.settings
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*rheo[^ \t]* //' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^[ \t]*include.*rheo.*$/d' ../Makefile.package.settings
fi
fi

10
src/RHEO/README Normal file
View File

@ -0,0 +1,10 @@
RHEO or Reproducing Hydrodynamics and Elastic Objects is a package to model
multiphase fluid systems. The authors include Joel Clemmer (Sandia), Thomas
O'Connor (Carnegie Mellon), and Eric Palermo (Carnegie Mellon).
Bond style rheo/shell, compute style rheo/property/atom, and fix style
rheo/temperature all depend on the BPM package.
This package requires the GNU scientific library (GSL). We recommend version
2.7 or later. To build this package, one must first separately install GSL in
a location that can be found by your environment.

174
src/RHEO/atom_vec_rheo.cpp Normal file
View File

@ -0,0 +1,174 @@
// clang-format off
/* ----------------------------------------------------------------------
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 authors:
Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "atom_vec_rheo.h"
#include "atom.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
AtomVecRHEO::AtomVecRHEO(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = Atom::ATOMIC;
mass_type = PER_TYPE;
forceclearflag = 1;
atom->rheo_status_flag = 1;
atom->pressure_flag = 1;
atom->rho_flag = 1;
atom->viscosity_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = {"rheo_status", "rho", "drho", "pressure", "viscosity"};
fields_copy = {"rheo_status", "rho", "drho", "pressure", "viscosity"};
fields_comm = {"rheo_status", "rho"};
fields_comm_vel = {"rheo_status", "rho"};
fields_reverse = {"drho"};
fields_border = {"rheo_status", "rho"};
fields_border_vel = {"rheo_status", "rho"};
fields_exchange = {"rheo_status", "rho"};
fields_restart = {"rheo_status", "rho"};
fields_create = {"rheo_status", "rho", "drho", "pressure", "viscosity"};
fields_data_atom = {"id", "type", "rheo_status", "rho", "x"};
fields_data_vel = {"id", "v"};
setup_fields();
}
/* ----------------------------------------------------------------------
set local copies of all grow ptrs used by this class, except defaults
needed in replicate when 2 atom classes exist and it calls pack_restart()
------------------------------------------------------------------------- */
void AtomVecRHEO::grow_pointers()
{
rheo_status = atom->rheo_status;
pressure = atom->pressure;
rho = atom->rho;
drho = atom->drho;
viscosity = atom->viscosity;
}
/* ----------------------------------------------------------------------
clear extra forces starting at atom N
nbytes = # of bytes to clear for a per-atom vector
------------------------------------------------------------------------- */
void AtomVecRHEO::force_clear(int n, size_t nbytes)
{
memset(&drho[n], 0, nbytes);
}
/* ----------------------------------------------------------------------
initialize non-zero atom quantities
------------------------------------------------------------------------- */
void AtomVecRHEO::create_atom_post(int ilocal)
{
rho[ilocal] = 1.0;
}
/* ----------------------------------------------------------------------
modify what AtomVec::data_atom() just unpacked
or initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecRHEO::data_atom_post(int ilocal)
{
drho[ilocal] = 0.0;
pressure[ilocal] = 0.0;
viscosity[ilocal] = 0.0;
}
/* ----------------------------------------------------------------------
assign an index to named atom property and return index
return -1 if name is unknown to this atom style
------------------------------------------------------------------------- */
int AtomVecRHEO::property_atom(const std::string &name)
{
if (name == "rheo_status") return 0;
if (name == "pressure") return 1;
if (name == "rho") return 2;
if (name == "drho") return 3;
if (name == "viscosity") return 4;
return -1;
}
/* ----------------------------------------------------------------------
pack per-atom data into buf for ComputePropertyAtom
index maps to data specific to this atom style
------------------------------------------------------------------------- */
void AtomVecRHEO::pack_property_atom(int index, double *buf, int nvalues, int groupbit)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
int n = 0;
if (index == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rheo_status[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = pressure[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 3) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = drho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 4) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = viscosity[i];
else
buf[n] = 0.0;
n += nvalues;
}
}
}

46
src/RHEO/atom_vec_rheo.h Normal file
View File

@ -0,0 +1,46 @@
/* -*- 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 ATOM_CLASS
// clang-format off
AtomStyle(rheo,AtomVecRHEO);
// clang-format on
#else
#ifndef LMP_ATOM_VEC_RHEO_H
#define LMP_ATOM_VEC_RHEO_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecRHEO : virtual public AtomVec {
public:
AtomVecRHEO(class LAMMPS *);
void grow_pointers() override;
void force_clear(int, size_t) override;
void create_atom_post(int) override;
void data_atom_post(int) override;
int property_atom(const std::string &) override;
void pack_property_atom(int, double *, int, int) override;
private:
int *rheo_status;
double *pressure, *rho, *drho, *viscosity;
};
} // namespace LAMMPS_NS
#endif
#endif

222
src/RHEO/atom_vec_rheo_thermal.cpp Normal file
View File

@ -0,0 +1,222 @@
// clang-format off
/* ----------------------------------------------------------------------
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 authors:
Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "atom_vec_rheo_thermal.h"
#include "atom.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = Atom::ATOMIC;
mass_type = PER_TYPE;
forceclearflag = 1;
atom->rheo_status_flag = 1;
atom->conductivity_flag = 1;
atom->temperature_flag = 1;
atom->esph_flag = 1;
atom->heatflow_flag = 1;
atom->pressure_flag = 1;
atom->rho_flag = 1;
atom->viscosity_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
fields_copy = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
fields_comm = {"rheo_status", "rho", "esph"};
fields_comm_vel = {"rheo_status", "rho", "esph"};
fields_reverse = {"drho", "heatflow"};
fields_border = {"rheo_status", "rho", "esph"};
fields_border_vel = {"rheo_status", "rho", "esph"};
fields_exchange = {"rheo_status", "rho", "esph"};
fields_restart = {"rheo_status", "rho", "esph"};
fields_create = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"};
fields_data_atom = {"id", "type", "rheo_status", "rho", "esph", "x"};
fields_data_vel = {"id", "v"};
setup_fields();
}
/* ----------------------------------------------------------------------
set local copies of all grow ptrs used by this class, except defaults
needed in replicate when 2 atom classes exist and it calls pack_restart()
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::grow_pointers()
{
rheo_status = atom->rheo_status;
conductivity = atom->conductivity;
temperature = atom->temperature;
esph = atom->esph;
heatflow = atom->heatflow;
pressure = atom->pressure;
rho = atom->rho;
drho = atom->drho;
viscosity = atom->viscosity;
}
/* ----------------------------------------------------------------------
clear extra forces starting at atom N
nbytes = # of bytes to clear for a per-atom vector
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::force_clear(int n, size_t nbytes)
{
memset(&drho[n], 0, nbytes);
memset(&heatflow[n], 0, nbytes);
}
/* ----------------------------------------------------------------------
initialize non-zero atom quantities
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::create_atom_post(int ilocal)
{
rho[ilocal] = 1.0;
}
/* ----------------------------------------------------------------------
modify what AtomVec::data_atom() just unpacked
or initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::data_atom_post(int ilocal)
{
drho[ilocal] = 0.0;
heatflow[ilocal] = 0.0;
temperature[ilocal] = 0.0;
pressure[ilocal] = 0.0;
viscosity[ilocal] = 0.0;
conductivity[ilocal] = 0.0;
}
/* ----------------------------------------------------------------------
assign an index to named atom property and return index
return -1 if name is unknown to this atom style
------------------------------------------------------------------------- */
int AtomVecRHEOThermal::property_atom(const std::string &name)
{
if (name == "rheo_status") return 0;
if (name == "rho") return 1;
if (name == "drho") return 2;
if (name == "temperature") return 3;
if (name == "esph") return 4;
if (name == "heatflow") return 5;
if (name == "conductivity") return 6;
if (name == "pressure") return 7;
if (name == "viscosity") return 8;
return -1;
}
/* ----------------------------------------------------------------------
pack per-atom data into buf for ComputePropertyAtom
index maps to data specific to this atom style
------------------------------------------------------------------------- */
void AtomVecRHEOThermal::pack_property_atom(int index, double *buf, int nvalues, int groupbit)
{
int *mask = atom->mask;
int nlocal = atom->nlocal;
int n = 0;
if (index == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rheo_status[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = rho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = drho[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 3) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = temperature[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 4) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = esph[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 5) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = heatflow[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 6) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = conductivity[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 7) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = pressure[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 8) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = viscosity[i];
else
buf[n] = 0.0;
n += nvalues;
}
}
}

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