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

Merge branch 'clean-master2' of github.com:julient31/lammps into gneb_spin

Browse Source
This commit is contained in:
julient31
2019-04-01 16:12:11 -06:00
parent a0e5d14990 585403d650
commit 677974573a
140 changed files with 11805 additions and 6184 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
doc/src/Build_extras.txt
View File

@ -82,17 +82,19 @@ which GPU hardware to build for.
[CMake build]:
-D GPU_API=value # value = opencl (default) or cuda
-D GPU_PREC=value # precision setting
# value = double or mixed (default) or single
-D OCL_TUNE=value # hardware choice for GPU_API=opencl
# generic (default) or intel (Intel CPU) or fermi, kepler, cypress (NVIDIA)
-D GPU_ARCH=value # primary GPU hardware choice for GPU_API=cuda
# value = sm_XX, see below
# default is Cuda-compiler dependent, but typically sm_20
-D CUDPP_OPT=value # optimization setting for GPU_API=cuda
# enables CUDA Performance Primitives Optimizations
# yes (default) or no :pre
-D GPU_API=value # value = opencl (default) or cuda
-D GPU_PREC=value # precision setting
# value = double or mixed (default) or single
-D OCL_TUNE=value # hardware choice for GPU_API=opencl
# generic (default) or intel (Intel CPU) or fermi, kepler, cypress (NVIDIA)
-D GPU_ARCH=value # primary GPU hardware choice for GPU_API=cuda
# value = sm_XX, see below
# default is Cuda-compiler dependent, but typically sm_20
-D CUDPP_OPT=value # optimization setting for GPU_API=cuda
# enables CUDA Performance Primitives Optimizations
# value = yes (default) or no
-D CUDA_MPS_SUPPORT=value # enables some tweaks required to run with active nvidia-cuda-mps daemon
# value = yes or no (default) :pre
GPU_ARCH settings for different GPU hardware is as follows:
@ -195,7 +197,7 @@ https://openkim.org/browse/models/by-model-drivers
If DOWNLOAD_KIM is set, the KIM library will be downloaded and built
inside the CMake build directory. If the KIM library is already on
your system (in a location CMake cannot find it), set the PKG_CONFIG_PATH
environment variable so that libkim-api-v2 can be found.
environment variable so that libkim-api can be found.
[Traditional make]:

2
doc/src/Commands_fix.txt
View File

@ -224,7 +224,7 @@ OPT.
"wall/body/polyhedron"_fix_wall_body_polyhedron.html,
"wall/colloid"_fix_wall.html,
"wall/ees"_fix_wall_ees.html,
"wall/gran (o)"_fix_wall_gran.html,
"wall/gran"_fix_wall_gran.html,
"wall/gran/region"_fix_wall_gran_region.html,
"wall/harmonic"_fix_wall.html,
"wall/lj1043"_fix_wall.html,

31
doc/src/Howto_pylammps.txt
View File

@ -57,6 +57,17 @@ library is then loaded by the Python interface. In this example we enable the
MOLECULE package and compile LAMMPS with C++ exceptions, PNG, JPEG and FFMPEG
output support enabled.
Step 1a: For the CMake based build system, the steps are:
mkdir $LAMMPS_DIR/build-shared
cd $LAMMPS_DIR/build-shared :pre
# MPI, PNG, Jpeg, FFMPEG are auto-detected
cmake ../cmake -DPKG_MOLECULE=yes -DLAMMPS_EXCEPTIONS=yes -DBUILD_LIB=yes -DBUILD_SHARED_LIBS=yes
make :pre
Step 1b: For the legacy, make based build system, the steps are:
cd $LAMMPS_DIR/src :pre
# add packages if necessary
@ -68,10 +79,9 @@ make mpi mode=shlib LMP_INC="-DLAMMPS_PNG -DLAMMPS_JPEG -DLAMMPS_FFMPEG -DLAMMPS
Step 2: Installing the LAMMPS Python package :h6
PyLammps is part of the lammps Python package. To install it simply install
that package into your current Python installation.
that package into your current Python installation with:
cd $LAMMPS_DIR/python
python install.py :pre
make install-python :pre
NOTE: Recompiling the shared library requires re-installing the Python package
@ -94,14 +104,21 @@ apt-get install python-virtualenv :pre
Creating a virtualenv with lammps installed :h6
# create virtualenv name 'testing' :pre
# create virtualenv named 'testing'
virtualenv $HOME/python/testing :pre
# activate 'testing' environment
source testing/bin/activate :pre
source $HOME/python/testing/bin/activate :pre
Now configure and compile the LAMMPS shared library as outlined above.
When using CMake and the shared library has already been build, you
need to re-run CMake to update the location of the python executable
to the location in the virtual environment with:
cmake . -DPYTHON_EXECUTABLE=$(which python) :pre
# install LAMMPS package in virtualenv
(testing) cd $LAMMPS_DIR/python
(testing) python install.py :pre
(testing) make install-python :pre
# install other useful packages
(testing) pip install matplotlib jupyter mpi4py :pre

4
doc/src/Manual.txt
View File

@ -1,7 +1,7 @@
<!-- HTML_ONLY -->
<HEAD>
<TITLE>LAMMPS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="28 Feb 2019 version">
<META NAME="docnumber" CONTENT="29 Mar 2019 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD>
@ -21,7 +21,7 @@
:line
LAMMPS Documentation :c,h1
28 Feb 2019 version :c,h2
29 Mar 2019 version :c,h2
"What is a LAMMPS version?"_Manual_version.html

69
doc/src/Python_install.txt
View File

@ -12,16 +12,23 @@ Installing LAMMPS in Python :h3
For Python to invoke LAMMPS, there are 2 files it needs to know about:
python/lammps.py
src/liblammps.so :ul
liblammps.so or liblammps.dylib :ul
Lammps.py is the Python wrapper on the LAMMPS library interface.
Liblammps.so is the shared LAMMPS library that Python loads, as
described above.
The python source code in lammps.py is the Python wrapper on the
LAMMPS library interface. The liblammps.so or liblammps.dylib file
is the shared LAMMPS library that Python loads dynamically.
You can insure Python can find these files in one of two ways:
You can achieve that Python can find these files in one of two ways:
set two environment variables
run the python/install.py script :ul
set two environment variables pointing to the location in the source tree
run "make install-python" or run the python/install.py script explicitly :ul
When calling "make install-python" LAMMPS will try to install the
python module and the shared library into the python site-packages folders;
either the system-wide ones, or the local users ones (in case of insufficient
permissions for the global install). Python will then find the module
and shared library file automatically. The exact location of these folders
depends on your python version and your operating system.
If you set the paths to these files as environment variables, you only
have to do it once. For the csh or tcsh shells, add something like
@ -30,42 +37,28 @@ this to your ~/.cshrc file, one line for each of the two files:
setenv PYTHONPATH $\{PYTHONPATH\}:/home/sjplimp/lammps/python
setenv LD_LIBRARY_PATH $\{LD_LIBRARY_PATH\}:/home/sjplimp/lammps/src :pre
If you use the python/install.py script, you need to invoke it every
time you rebuild LAMMPS (as a shared library) or make changes to the
python/lammps.py file.
On MacOSX you may also need to set DYLD_LIBRARY_PATH accordingly.
For Bourne/Korn shells accordingly into the corresponding files using
the "export" shell builtin.
You can invoke install.py from the python directory as
If you use "make install-python" or the python/install.py script, you need
to invoke it every time you rebuild LAMMPS (as a shared library) or
make changes to the python/lammps.py file, so that the site-packages
files are updated with the new version.
% python install.py \[libdir\] \[pydir\] :pre
If the default settings of "make install-python" are not what you want,
you can invoke install.py from the python directory manually as
The optional libdir is where to copy the LAMMPS shared library to; the
default is /usr/local/lib. The optional pydir is where to copy the
lammps.py file to; the default is the site-packages directory of the
version of Python that is running the install script.
% python install.py -m \<python module\> -l <shared library> -v <version.h file> \[-d \<pydir\>\] :pre
Note that libdir must be a location that is in your default
LD_LIBRARY_PATH, like /usr/local/lib or /usr/lib. And pydir must be a
location that Python looks in by default for imported modules, like
its site-packages dir. If you want to copy these files to
non-standard locations, such as within your own user space, you will
need to set your PYTHONPATH and LD_LIBRARY_PATH environment variables
accordingly, as above.
The -m flag points to the lammps.py python module file to be installed,
the -l flag points to the LAMMPS shared library file to be installed,
the -v flag points to the version.h file in the LAMMPS source
and the optional -d flag to a custom (legacy) installation folder :ul
If the install.py script does not allow you to copy files into system
directories, prefix the python command with "sudo". If you do this,
make sure that the Python that root runs is the same as the Python you
run. E.g. you may need to do something like
% sudo /usr/local/bin/python install.py \[libdir\] \[pydir\] :pre
You can also invoke install.py from the make command in the src
directory as
% make install-python :pre
In this mode you cannot append optional arguments. Again, you may
need to prefix this with "sudo". In this mode you cannot control
which Python is invoked by root.
If you use a legacy installation folder, you will need to set your
PYTHONPATH and LD_LIBRARY_PATH (and/or DYLD_LIBRARY_PATH) environment
variables accordingly, as described above.
Note that if you want Python to be able to load different versions of
the LAMMPS shared library (see "this section"_Python_shlib.html), you will

19
doc/src/Python_overview.txt
View File

@ -13,11 +13,11 @@ Overview of Python and LAMMPS :h3
LAMMPS can work together with Python in three ways. First, Python can
wrap LAMMPS through the its "library interface"_Howto_library.html, so
that a Python script can create one or more instances of LAMMPS and
launch one or more simulations. In Python lingo, this is "extending"
Python with LAMMPS.
launch one or more simulations. In Python lingo, this is called
"extending" Python with a LAMMPS module.
Second, a lower-level Python interface can be used indirectly through
provided PyLammps and IPyLammps wrapper classes, written in Python.
the provided PyLammps and IPyLammps wrapper classes, written in Python.
These wrappers try to simplify the usage of LAMMPS in Python by
providing an object-based interface to common LAMMPS functionality.
They also reduces the amount of code necessary to parameterize LAMMPS
@ -25,11 +25,12 @@ scripts through Python and make variables and computes directly
accessible.
Third, LAMMPS can use the Python interpreter, so that a LAMMPS
input script can invoke Python code directly, and pass information
back-and-forth between the input script and Python functions you
write. This Python code can also callback to LAMMPS to query or change
its attributes. In Python lingo, this is "embedding" Python in
LAMMPS. When used in this mode, Python can perform operations that
the simple LAMMPS input script syntax cannot.
input script or styles can invoke Python code directly, and pass
information back-and-forth between the input script and Python
functions you write. This Python code can also callback to LAMMPS
to query or change its attributes through the LAMMPS Python module
mentioned above. In Python lingo, this is "embedding" Python in
LAMMPS. When used in this mode, Python can perform script operations
that the simple LAMMPS input script syntax can not.

4
doc/src/Speed_intel.txt
View File

@ -34,6 +34,10 @@ rebo, sw, tersoff :l
K-Space Styles: pppm, pppm/disp :l
:ule
IMPORTANT NOTE: None of the styles in the USER-INTEL package currently
support computing per-atom stress. If any compute or fix in your
input requires it, LAMMPS will abort with an error message.
[Speed-ups to expect:]
The speedups will depend on your simulation, the hardware, which

7
doc/src/compute.txt
View File

@ -54,9 +54,10 @@ local quantities have the word "local" in their style,
e.g. {bond/local}. Styles with neither "atom" or "local" in their
style produce global quantities.
Note that a single compute produces either global or per-atom or local
quantities, but never more than one of these (with only a few
exceptions, as documented by individual compute commands).
Note that a single compute can produce either global or per-atom or
local quantities, but not both global and per-atom. It can produce
local quantities in tandem with global or per-atom quantities. The
compute doc page will explain.
Global, per-atom, and local quantities each come in three kinds: a
single scalar value, a vector of values, or a 2d array of values. The

6
doc/src/fix.txt
View File

@ -83,8 +83,10 @@ not in the specified fix group. Local quantities are calculated by
each processor based on the atoms it owns, but there may be zero or
more per atoms.
Note that a single fix may produces either global or per-atom or local
quantities (or none at all), but never more than one of these.
Note that a single fix can produce either global or per-atom or local
quantities (or none at all), but not both global and per-atom. It can
produce local quantities in tandem with global or per-atom quantities.
The fix doc page will explain.
Global, per-atom, and local quantities each come in three kinds: a
single scalar value, a vector of values, or a 2d array of values. The

19
doc/src/fix_ave_histo.txt
View File

@ -35,6 +35,7 @@ keyword = {mode} or {file} or {ave} or {start} or {beyond} or {overwrite} or {ti
{mode} arg = {scalar} or {vector}
scalar = all input values are scalars
vector = all input values are vectors
{kind} arg = {global} or {peratom} or {local}
{file} arg = filename
filename = name of file to output histogram(s) to
{ave} args = {one} or {running} or {window}
@ -92,7 +93,8 @@ either all global, all per-atom, or all local quantities. Inputs of
different kinds (e.g. global and per-atom) cannot be mixed. Atom
attributes are per-atom vector values. See the doc page for
individual "compute" and "fix" commands to see what kinds of
quantities they generate.
quantities they generate. See the optional {kind} keyword below for
how to force the fix ave/histo command to disambiguate if necessary.
Note that the output of this command is a single histogram for all
input values combined together, not one histogram per input value.
@ -231,6 +233,14 @@ keyword is set to {vector}, then all input values must be global or
per-atom or local vectors, or columns of global or per-atom or local
arrays.
The {kind} keyword only needs to be set if a compute or fix produces
more than one kind of output (global, per-atom, local). If this is
not the case, then LAMMPS will determine what kind of input is
provided and whether all the input arguments are consistent. If a
compute or fix produces more than one kind of output, the {kind}
keyword should be used to specify which output will be used. The
remaining input arguments must still be consistent.
The {beyond} keyword determines how input values that fall outside the
{lo} to {hi} bounds are treated. Values such that {lo} <= value <=
{hi} are assigned to one bin. Values on a bin boundary are assigned
@ -240,7 +250,7 @@ If {beyond} is set to {end} then values < {lo} are counted in the
first bin and values > {hi} are counted in the last bin. If {beyond}
is set to {extend} then two extra bins are created, so that there are
Nbins+2 total bins. Values < {lo} are counted in the first bin and
values > {hi} are counted in the last bin (Nbins+1). Values between
values > {hi} are counted in the last bin (Nbins+2). Values between
{lo} and {hi} (inclusive) are counted in bins 2 through Nbins+1. The
"coordinate" stored and printed for these two extra bins is {lo} and
{hi}.
@ -354,5 +364,6 @@ ave/chunk"_fix_ave_chunk.html, "fix ave/time"_fix_ave_time.html,
[Default:] none
The option defaults are mode = scalar, ave = one, start = 0, no file
output, beyond = ignore, and title 1,2,3 = strings as described above.
The option defaults are mode = scalar, kind = figured out from input
arguments, ave = one, start = 0, no file output, beyond = ignore, and
title 1,2,3 = strings as described above.

45
doc/src/fix_hyper_global.txt
View File

@ -102,7 +102,7 @@ Bi = exp(beta * Vij(max)) :pre
where beta = 1/kTequil, and {Tequil} is the temperature of the system
and an argument to this fix. Note that Bi >= 1 at every step.
NOTE: To run GHD, the input script must also use the "fix
NOTE: To run a GHD simulation, the input script must also use the "fix
langevin"_fix_langevin.html command to thermostat the atoms at the
same {Tequil} as specified by this fix, so that the system is running
constant-temperature (NVT) dynamics. LAMMPS does not check that this
@ -166,9 +166,9 @@ correctly. There will just be fewer events because the hyper time
NOTE: If you have no physical intuition as to the smallest barrier
height in your system, a reasonable strategy to determine the largest
{Vmax} you can use for an LHD model, is to run a sequence of
{Vmax} you can use for a GHD model, is to run a sequence of
simulations with smaller and smaller {Vmax} values, until the event
rate does not change.
rate does not change (as a function of hyper time).
The {Tequil} argument is the temperature at which the system is
simulated; see the comment above about the "fix
@ -177,7 +177,8 @@ beta term in the exponential factor that determines how much boost is
achieved as a function of the bias potential.
In general, the lower the value of {Tequil} and the higher the value
of {Vmax}, the more boost will be achievable by the GHD algorithm.
of {Vmax}, the more time boost will be achievable by the GHD
algorithm.
:line
@ -190,41 +191,43 @@ The "fix_modify"_fix_modify.html {energy} option is supported by this
fix to add the energy of the bias potential to the the system's
potential energy as part of "thermodynamic output"_thermo_style.html.
This fix computes a global scalar and global vector of length 11, which
This fix computes a global scalar and global vector of length 12, which
can be accessed by various "output commands"_Howto_output.html. The
scalar is the magnitude of the bias potential (energy units) applied on
the current timestep. The vector stores the following quantities:
1 = boost factor on this step (unitless)
2 = max strain Eij of any bond on this step (unitless)
2 = max strain Eij of any bond on this step (absolute value, unitless)
3 = ID of first atom in the max-strain bond
4 = ID of second atom in the max-strain bond
5 = average # of bonds/atom on this step :ul
6 = fraction of timesteps with bias = 0.0 during this run
7 = max drift distance of any atom during this run (distance units)
8 = max bond length during this run (distance units) :ul
6 = fraction of timesteps where the biased bond has bias = 0.0 during this run
7 = fraction of timesteps where the biased bond has negative strain during this run
8 = max drift distance of any atom during this run (distance units)
9 = max bond length during this run (distance units) :ul
9 = cumulative hyper time since fix was defined (time units)
10 = cumulative count of event timesteps since fix was defined
11 = cumulative count of atoms in events since fix was defined :ul
10 = cumulative hyper time since fix was defined (time units)
11 = cumulative count of event timesteps since fix was defined
12 = cumulative count of atoms in events since fix was defined :ul
The first 5 quantities are for the current timestep. Quantities 6-8
are for the current hyper run. Quantities 9-11 are cumulative across
multiple runs (since the fix was defined in the input script).
The first 5 quantities are for the current timestep. Quantities 6-9
are for the current hyper run. They are reset each time a new hyper
run is performed. Quantities 19-12 are cumulative across multiple
runs (since the point in the input script the fix was defined).
For value 7, drift is the distance an atom moves between timesteps
when the bond list is reset, i.e. between events. Atoms involved in
an event will typically move the greatest distance since others are
typically oscillating around their lattice site.
For value 8, drift is the distance an atom moves between two quenched
states when the second quench determines an event has occurred. Atoms
involved in an event will typically move the greatest distance since
others typically remain near their original quenched position.
For value 10, events are checked for by the "hyper"_hyper.html command
For value 11, events are checked for by the "hyper"_hyper.html command
once every {Nevent} timesteps. This value is the count of those
timesteps on which one (or more) events was detected. It is NOT the
number of distinct events, since more than one event may occur in the
same {Nevent} time window.
For value 11, each time the "hyper"_hyper.html command checks for an
For value 12, each time the "hyper"_hyper.html command checks for an
event, it invokes a compute to flag zero or more atoms as
participating in one or more events. E.g. atoms that have displaced
more than some distance from the previous quench state. Value 11 is

225
doc/src/fix_hyper_local.txt
View File

@ -22,10 +22,9 @@ Dcut = minimum distance between boosted bonds (distance units) :l
alpha = boostostat relaxation time (time units) :l
Btarget = desired time boost factor (unitless) :l
zero or more keyword/value pairs may be appended :l
keyword = {lost} or {check/bias} or {check/coeff}
{lostbond} value = error/warn/ignore
{check/bias} values = Nevery error/warn/ignore
{check/coeff} values = Nevery error/warn/ignore :pre
keyword = {check/ghost} or {check/bias} :l
{check/ghost} values = none
{check/bias} values = Nevery error/warn/ignore :pre
:ule
[Examples:]
@ -65,8 +64,8 @@ To understand this description, you should first read the description
of the GHD algorithm on the "fix hyper/global"_fix_hyper_global.html
doc page. This description of LHD builds on the GHD description.
The definition of bonds, Eij, and Emax are the same for GHD and LHD.
The formulas for Vij(max) and Fij(max) are also the same except for a
The definition of bonds and Eij are the same for GHD and LHD. The
formulas for Vij(max) and Fij(max) are also the same except for a
pre-factor Cij, explained below.
The bias energy Vij applied to a bond IJ with maximum strain is
@ -117,11 +116,11 @@ where Vkl(max) is the bias energy of the maxstrain bond KL within bond
IJ's neighborhood, beta = 1/kTequil, and {Tequil} is the temperature
of the system and an argument to this fix.
NOTE: To run LHD, the input script must also use the "fix
langevin"_fix_langevin.html command to thermostat the atoms at the
same {Tequil} as specified by this fix, so that the system is running
constant-temperature (NVT) dynamics. LAMMPS does not check that this
is done.
NOTE: To run an LHD simulation, the input script must also use the
"fix langevin"_fix_langevin.html command to thermostat the atoms at
the same {Tequil} as specified by this fix, so that the system is
running constant-temperature (NVT) dynamics. LAMMPS does not check
that this is done.
Note that if IJ = KL, then bond IJ is a biased bond on that timestep,
otherwise it is not. But regardless, the boost factor Bij can be
@ -216,20 +215,20 @@ each pair. E.g. something like 2x the cutoff of the interatomic
potential. In practice a {Dcut} value of ~10 Angstroms seems to work
well for many solid-state systems.
NOTE: You must also insure that ghost atom communication is performed
for a distance of at least {Dcut} + {cutevent} where {cutevent} = the
distance one or more atoms move (between quenched states) to be
considered an "event". It is an argument to the "compute
event/displace" command used to detect events. By default the ghost
communication distance is set by the pair_style cutoff, which will
typically be < {Dcut}. The "comm_modify cutoff"_comm_modify.html
command can be used to set the ghost cutoff explicitly, e.g.
NOTE: You should insure that ghost atom communication is performed for
a distance of at least {Dcut} + {cutevent} = the distance one or more
atoms move (between quenched states) to be considered an "event". It
is an argument to the "compute event/displace" command used to detect
events. By default the ghost communication distance is set by the
pair_style cutoff, which will typically be < {Dcut}. The "comm_modify
cutoff"_comm_modify.html command should be used to override the ghost
cutoff explicitly, e.g.
comm_modify cutoff 12.0 :pre
This fix does not know the {cutevent} parameter, but uses half the
bond length as an estimate to warn if the ghost cutoff is not long
enough.
Note that this fix does not know the {cutevent} parameter, but uses
half the {cutbond} parameter as an estimate to warn if the ghost
cutoff is not long enough.
As described above the {alpha} argument is a pre-factor in the
boostostat update equation for each bond's Cij prefactor. {Alpha} is
@ -269,7 +268,30 @@ NOTE: If you have no physical intuition as to the smallest barrier
height in your system, a reasonable strategy to determine the largest
{Btarget} you can use for an LHD model, is to run a sequence of
simulations with smaller and smaller {Btarget} values, until the event
rate does not change.
rate does not change (as a function of hyper time).
:line
Here is additional information on the optional keywords for this fix.
The {check/ghost} keyword turns on extra computation each timestep to
compute statistics about ghost atoms used to determine which bonds to
bias. The output of these stats are the vector values 14 and 15,
described below. If this keyword is not enabled, the output
of the stats will be zero.
The {check/bias} keyword turns on extra computation and communication
to check if any biased bonds are closer than {Dcut} to each other,
which should not be the case if LHD is operating correctly. Thus it
is a debugging check. The {Nevery} setting determines how often the
check is made. The {error}, {warn}, or {ignore} setting determines
what is done if the count of too-close bonds is not zero. Either the
code will exit, or issue a warning, or silently tally the count. The
count can be output as vector value 17, as described below. If this
keyword is not enabled, the output of that statistic will be 0.
Note that both of these computations are costly, hence they are only
enabled by these keywords.
:line
@ -282,95 +304,120 @@ The "fix_modify"_fix_modify.html {energy} option is supported by this
fix to add the energy of the bias potential to the the system's
potential energy as part of "thermodynamic output"_thermo_style.html.
This fix computes a global scalar and global vector of length 23,
which can be accessed by various "output
commands"_Howto_output.html. The scalar is the magnitude of
the bias potential (energy units) applied on the current timestep,
summed over all biased bonds. The vector stores the following
quantities:
This fix computes a global scalar and global vector of length 21,
which can be accessed by various "output commands"_Howto_output.html.
The scalar is the magnitude of the bias potential (energy units)
applied on the current timestep, summed over all biased bonds. The
vector stores the following quantities:
1 = # of biased bonds on this step
2 = max strain Eij of any bond on this step (unitless)
3 = average bias potential for all biased bonds on this step (energy units)
2 = max strain Eij of any bond on this step (absolute value, unitless)
3 = average bias coeff for all bonds on this step (unitless)
4 = average # of bonds/atom on this step
5 = average neighbor bonds/bond on this step within {Dcut} :ul
6 = fraction of steps and bonds with no bias during this run
7 = max drift distance of any atom during this run (distance units)
8 = max bond length during this run (distance units)
9 = average # of biased bonds/step during this run
10 = average bias potential for all biased bonds during this run (energy units)
11 = max bias potential for any biased bond during this run (energy units)
12 = min bias potential for any biased bond during this run (energy units)
13 = max distance from my sub-box of any ghost atom with maxstrain < qfactor during this run (distance units)
14 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units)
15 = count of ghost neighbor atoms not found on reneighbor steps during this run
16 = count of lost bond partners during this run
17 = average bias coeff for lost bond partners during this run
18 = count of bias overlaps found during this run
19 = count of non-matching bias coefficients found during this run :ul
6 = max bond length during this run (distance units)
7 = average # of biased bonds/step during this run
8 = fraction of biased bonds with no bias during this run
9 = fraction of biased bonds with negative strain during this run
10 = average bias coeff for all bonds during this run (unitless)
11 = min bias coeff for any bond during this run (unitless)
12 = max bias coeff for any bond during this run (unitless)
20 = cumulative hyper time since fix created (time units)
21 = cumulative count of event timesteps since fix created
22 = cumulative count of atoms in events since fix created
23 = cumulative # of new bonds since fix created :ul
13 = max drift distance of any bond atom during this run (distance units)
14 = max distance from proc subbox of any ghost atom with maxstrain < qfactor during this run (distance units)
15 = max distance outside my box of any ghost atom with any maxstrain during this run (distance units)
16 = count of ghost atoms that could not be found on reneighbor steps during this run
17 = count of bias overlaps (< Dcut) found during this run
18 = cumulative hyper time since fix created (time units)
19 = cumulative count of event timesteps since fix created
20 = cumulative count of atoms in events since fix created
21 = cumulative # of new bonds formed since fix created :ul
The first quantities (1-5) are for the current timestep. Quantities
6-19 are for the current hyper run. They are reset each time a new
hyper run is performed. Quantities 20-23 are cumulative across
multiple runs (since the fix was defined in the input script).
6-17 are for the current hyper run. They are reset each time a new
hyper run is performed. Quantities 18-21 are cumulative across
multiple runs (since the point in the input script the fix was
defined).
For value 6, the numerator is a count of all biased bonds on every
For value 8, the numerator is a count of all biased bonds on each
timestep whose bias energy = 0.0 due to Eij >= {qfactor}. The
denominator is the count of all biased bonds on all timesteps.
For value 7, drift is the distance an atom moves between timesteps
when the bond list is reset, i.e. between events. Atoms involved in
an event will typically move the greatest distance since others are
typically oscillating around their lattice site.
For value 9, the numerator is a count of all biased bonds on each
timestep with negative strain. The denominator is the count of all
biased bonds on all timesteps.
For values 13 and 14, the maxstrain of a ghost atom is the maxstrain
of any bond it is part of, and it is checked for ghost atoms within
the bond neighbor cutoff.
Values 13-17 are mostly useful for debugging and diagnostic purposes.
Values 15-19 are mostly useful for debugging and diagnostic purposes.
For value 13, drift is the distance an atom moves between two quenched
states when the second quench determines an event has occurred. Atoms
involved in an event will typically move the greatest distance since
others typically remain near their original quenched position.
For values 15-17, it is possible that a ghost atom owned by another
processor will move far enough (e.g. as part of an event-in-progress)
that it will no longer be within the communication cutoff distance for
acquiring ghost atoms. Likewise it may be a ghost atom bond partner
that cannot be found because it has moved too far. These values count
those occurrences. Because they typically involve atoms that are part
of events, they do not usually indicate bad dynamics. Value 16 is the
average bias coefficient for bonds where a partner atom was lost.
For values 14-16, neighbor atoms in the full neighbor list with cutoff
{Dcut} may be ghost atoms outside a processor's sub-box. Before the
next event occurs they may move further than {Dcut} away from the
sub-box boundary. Value 14 is the furthest (from the sub-box) any
ghost atom in the neighbor list with maxstrain < {qfactor} was
accessed during the run. Value 15 is the same except that the ghost
atom's maxstrain may be >= {qfactor}, which may mean it is about to
participate in an event. Value 16 is a count of how many ghost atoms
could not be found on reneighbor steps, presumably because they moved
too far away due to their participation in an event (which will likely
be detected at the next quench).
For value 18, no two bonds should be biased if they are within a
Typical values for 14 and 15 should be slightly larger than {Dcut},
which accounts for ghost atoms initially at a {Dcut} distance moving
thermally before the next event takes place.
Note that for values 14 and 15 to be computed, the optional keyword
{check/ghost} must be specified. Otherwise these values will be zero.
This is because computing them incurs overhead, so the values are only
computed if requested.
Value 16 should be zero or small. As explained above a small count
likely means some ghost atoms were participating in their own events
and moved a longer distance. If the value is large, it likely means
the communication cutoff for ghosts is too close to {Dcut} leading to
many not-found ghost atoms before the next event. This may lead to a
reduced number of bonds being selected for biasing, since the code
assumes those atoms are part of highly strained bonds. As explained
above, the "comm_modify cutoff"_comm_modify.html command can be used
to set a longer cutoff.
For value 17, no two bonds should be biased if they are within a
{Dcut} distance of each other. This value should be zero, indicating
that no pair of bonds "overlap", meaning they are closer than {Dcut}
from each other.
that no pair of biased bonds are closer than {Dcut} from each other.
For value 19, the same bias coefficient is stored by both atoms in an
IJ bond. This value should be zero, indicating that for all bonds,
each atom in the bond stores the a bias coefficient with the same
value.
Note that for values 17 to be computed, the optional keyword
{check/bias} must be specified and it determines how often this check
is performed. This is because performing the check incurs overhead,
so if only computed as often as requested.
Value 20 is simply the specified {boost} factor times the number of
timestep times the timestep size.
The result at the end of the run is the cumulative total from every
timestep the check was made. Note that the value is a count of atoms
in bonds which found other atoms in bonds too close, so it is almost
always an over-count of the number of too-close bonds.
For value 21, events are checked for by the "hyper"_hyper.html command
Value 18 is simply the specified {boost} factor times the number of
timesteps times the timestep size.
For value 19, events are checked for by the "hyper"_hyper.html command
once every {Nevent} timesteps. This value is the count of those
timesteps on which one (or more) events was detected. It is NOT the
number of distinct events, since more than one event may occur in the
same {Nevent} time window.
For value 22, each time the "hyper"_hyper.html command checks for an
For value 20, each time the "hyper"_hyper.html command checks for an
event, it invokes a compute to flag zero or more atoms as
participating in one or more events. E.g. atoms that have displaced
more than some distance from the previous quench state. Value 22 is
more than some distance from the previous quench state. Value 20 is
the cumulative count of the number of atoms participating in any of
the events that were found.
Value 23 tallies the number of new bonds created by the bond reset
Value 21 tallies the number of new bonds created by the bond reset
operation. Bonds between a specific I,J pair of atoms may persist for
the entire hyperdynamics simulation if neither I or J are involved in
an event.
@ -378,6 +425,16 @@ an event.
The scalar and vector values calculated by this fix are all
"intensive".
This fix also computes a local vector of length the number of bonds
currently in the system. The value for each bond is its Cij prefactor
(bias coefficient). These values can be can be accessed by various
"output commands"_Howto_output.html. A particularly useful one is the
"fix ave/histo"_fix_ave_histo.html command which can be used to
histogram the Cij values to see if they are distributed reasonably
close to 1.0, which indicates a good choice of {Vmax}.
The local values calculated by this fix are unitless.
No parameter of this fix can be used with the {start/stop} keywords of
the "run"_run.html command. This fix is not invoked during "energy
minimization"_minimize.html.
@ -392,7 +449,9 @@ doc page for more info.
"hyper"_hyper.html, "fix hyper/global"_fix_hyper_global.html
[Default:] None
[Default:]
The check/ghost and check/bias keywords are not enabled by default.
:line

101
doc/src/fix_wall_gran.txt
View File

@ -7,22 +7,24 @@
:line
fix wall/gran command :h3
fix wall/gran/omp command :h3
[Syntax:]
fix ID group-ID wall/gran fstyle Kn Kt gamma_n gamma_t xmu dampflag wallstyle args keyword values ... :pre
fix ID group-ID wall/gran fstyle fstyle_params wallstyle args keyword values ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
wall/gran = style name of this fix command :l
fstyle = style of force interactions between particles and wall :l
possible choices: hooke, hooke/history, hertz/history :pre
Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) :l
Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) :l
gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) :l
gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) :l
xmu = static yield criterion (unitless value between 0.0 and 1.0e4) :l
dampflag = 0 or 1 if tangential damping force is excluded or included :l
possible choices: hooke, hooke/history, hertz/history, granular :pre
fstyle_params = parameters associated with force interaction style :l
For {hooke}, {hooke/history}, and {hertz/history}, {fstyle_params} are:
Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below)
Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below)
gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below)
gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below)
xmu = static yield criterion (unitless value between 0.0 and 1.0e4)
dampflag = 0 or 1 if tangential damping force is excluded or included :pre
For {granular}, {fstyle_params} are set using the same syntax as for the {pair_coeff} command of "pair_style granular"_pair_granular.html :pre
wallstyle = {xplane} or {yplane} or {zplane} or {zcylinder} :l
args = list of arguments for a particular style :l
{xplane} or {yplane} or {zplane} args = lo hi
@ -44,7 +46,10 @@ keyword = {wiggle} or {shear} :l
fix 1 all wall/gran hooke 200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0
fix 1 all wall/gran hooke/history 200000.0 NULL 50.0 NULL 0.5 0 zplane 0.0 NULL
fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0 :pre
fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0
fix 3 all wall/gran granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 zplane 0.0 NULL
fix 4 all wall/gran granular jkr 1000.0 50.0 0.3 5.0 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall zcylinder 15.0 wiggle z 3.0 2.0
fix 5 all wall/gran granular dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall zplane 0.0 NULL :pre
[Description:]
@ -54,31 +59,40 @@ close enough to touch it.
The nature of the wall/particle interactions are determined by the
{fstyle} setting. It can be any of the styles defined by the
"pair_style granular"_pair_gran.html commands. Currently this is
{hooke}, {hooke/history}, or {hertz/history}. The equation for the
force between the wall and particles touching it is the same as the
corresponding equation on the "pair_style granular"_pair_gran.html doc
page, in the limit of one of the two particles going to infinite
radius and mass (flat wall). Specifically, delta = radius - r =
overlap of particle with wall, m_eff = mass of particle, and the
effective radius of contact = RiRj/Ri+Rj is just the radius of the
particle.
"pair_style gran/*"_pair_gran.html or the more general "pair_style
granular"_pair_granular.html" commands. Currently the options are
{hooke}, {hooke/history}, or {hertz/history} for the former, and
{granular} with all the possible options of the associated
{pair_coeff} command for the latter. The equation for the force
between the wall and particles touching it is the same as the
corresponding equation on the "pair_style gran/*"_pair_gran.html and
"pair_style_granular"_pair_granular.html doc pages, in the limit of
one of the two particles going to infinite radius and mass (flat
wall). Specifically, delta = radius - r = overlap of particle with
wall, m_eff = mass of particle, and the effective radius of contact =
RiRj/Ri+Rj is set to the radius of the particle.
The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag}
have the same meaning and units as those specified with the
"pair_style granular"_pair_gran.html commands. This means a NULL can
be used for either {Kt} or {gamma_t} as described on that page. If a
"pair_style gran/*"_pair_gran.html commands. This means a NULL can be
used for either {Kt} or {gamma_t} as described on that page. If a
NULL is used for {Kt}, then a default value is used where {Kt} = 2/7
{Kn}. If a NULL is used for {gamma_t}, then a default value is used
where {gamma_t} = 1/2 {gamma_n}.
All the model choices for cohesion, tangential friction, rolling
friction and twisting friction supported by the "pair_style
granular"_pair_granular.html through its {pair_coeff} command are also
supported for walls. These are discussed in greater detail on the doc
page for "pair_style granular"_pair_granular.html.
Note that you can choose a different force styles and/or different
values for the 6 wall/particle coefficients than for particle/particle
values for the wall/particle coefficients than for particle/particle
interactions. E.g. if you wish to model the wall as a different
material.
NOTE: As discussed on the doc page for "pair_style
granular"_pair_gran.html, versions of LAMMPS before 9Jan09 used a
gran/*"_pair_gran.html, versions of LAMMPS before 9Jan09 used a
different equation for Hertzian interactions. This means Hertizian
wall/particle interactions have also changed. They now include a
sqrt(radius) term which was not present before. Also the previous
@ -108,14 +122,14 @@ Optionally, the wall can be moving, if the {wiggle} or {shear}
keywords are appended. Both keywords cannot be used together.
For the {wiggle} keyword, the wall oscillates sinusoidally, similar to
the oscillations of particles which can be specified by the
"fix move"_fix_move.html command. This is useful in packing
simulations of granular particles. The arguments to the {wiggle}
keyword specify a dimension for the motion, as well as it's
{amplitude} and {period}. Note that if the dimension is in the plane
of the wall, this is effectively a shearing motion. If the dimension
is perpendicular to the wall, it is more of a shaking motion. A
{zcylinder} wall can only be wiggled in the z dimension.
the oscillations of particles which can be specified by the "fix
move"_fix_move.html command. This is useful in packing simulations of
granular particles. The arguments to the {wiggle} keyword specify a
dimension for the motion, as well as it's {amplitude} and {period}.
Note that if the dimension is in the plane of the wall, this is
effectively a shearing motion. If the dimension is perpendicular to
the wall, it is more of a shaking motion. A {zcylinder} wall can only
be wiggled in the z dimension.
Each timestep, the position of a wiggled wall in the appropriate {dim}
is set according to this equation:
@ -137,28 +151,6 @@ the clockwise direction for {vshear} > 0 or counter-clockwise for
{vshear} < 0. In this case, {vshear} is the tangential velocity of
the wall at whatever {radius} has been defined.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed on the "Speed packages"_Speed_packages.html doc
page. The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.
These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively. They are only enabled if
LAMMPS was built with those packages. See the "Build
package"_Build_package.html doc page for more info.
You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the "-suffix command-line
switch"_Run_options.html when you invoke LAMMPS, or you can use the
"suffix"_suffix.html command in your input script.
See the "Speed packages"_Speed_packages.html doc page for more
instructions on how to use the accelerated styles effectively.
[Restart, fix_modify, output, run start/stop, minimize info:]
This fix writes the shear friction state of atoms interacting with the
@ -188,6 +180,7 @@ Any dimension (xyz) that has a granular wall must be non-periodic.
"fix move"_fix_move.html,
"fix wall/gran/region"_fix_wall_gran_region.html,
"pair_style granular"_pair_gran.html
"pair_style gran/*"_pair_gran.html
"pair_style granular"_pair_granular.html
[Default:] none

76
doc/src/fix_wall_gran_region.txt
View File

@ -10,24 +10,30 @@ fix wall/gran/region command :h3
[Syntax:]
fix ID group-ID wall/gran/region fstyle Kn Kt gamma_n gamma_t xmu dampflag wallstyle regionID :pre
fix ID group-ID wall/gran/region fstyle fstyle_params wallstyle regionID :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
wall/region = style name of this fix command :l
fstyle = style of force interactions between particles and wall :l
possible choices: hooke, hooke/history, hertz/history :pre
Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below) :l
Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below) :l
gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below) :l
gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below) :l
xmu = static yield criterion (unitless value between 0.0 and 1.0e4) :l
dampflag = 0 or 1 if tangential damping force is excluded or included :l
possible choices: hooke, hooke/history, hertz/history, granular :pre
fstyle_params = parameters associated with force interaction style :l
For {hooke}, {hooke/history}, and {hertz/history}, {fstyle_params} are:
Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below)
Kt = elastic constant for tangential contact (force/distance units or pressure units - see discussion below)
gamma_n = damping coefficient for collisions in normal direction (1/time units or 1/time-distance units - see discussion below)
gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below)
xmu = static yield criterion (unitless value between 0.0 and 1.0e4)
dampflag = 0 or 1 if tangential damping force is excluded or included :pre
For {granular}, {fstyle_params} are set using the same syntax as for the {pair_coeff} command of "pair_style granular"_pair_granular.html :pre
wallstyle = region (see "fix wall/gran"_fix_wall_gran.html for options for other kinds of walls) :l
region-ID = region whose boundary will act as wall :l,ule
[Examples:]
fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone :pre
fix wall all wall/gran/region hooke/history 1000.0 200.0 200.0 100.0 0.5 1 region myCone
fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 region myBox
fix 4 all wall/gran/region granular jkr 1000.0 50.0 tangential linear_history 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone
fix 5 all wall/gran/region granular dmt 1000.0 50.0 0.3 10.0 tangential linear_history 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall region myCone :pre
[Description:]
@ -42,8 +48,8 @@ Here are snapshots of example models using this command.
Corresponding input scripts can be found in examples/granregion.
Click on the images to see a bigger picture. Movies of these
simulations are "here on the Movies
page"_http://lammps.sandia.gov/movies.html#granregion of the
LAMMPS web site.
page"_http://lammps.sandia.gov/movies.html#granregion of the LAMMPS
web site.
:image(JPG/gran_funnel_small.jpg,JPG/gran_funnel.png)
:image(JPG/gran_mixer_small.jpg,JPG/gran_mixer.png)
@ -123,12 +129,16 @@ to make the two faces differ by epsilon in their position.
The nature of the wall/particle interactions are determined by the
{fstyle} setting. It can be any of the styles defined by the
"pair_style granular"_pair_gran.html commands. Currently this is
{hooke}, {hooke/history}, or {hertz/history}. The equation for the
force between the wall and particles touching it is the same as the
corresponding equation on the "pair_style granular"_pair_gran.html doc
page, but the effective radius is calculated using the radius of the
particle and the radius of curvature of the wall at the contact point.
"pair_style gran/*"_pair_gran.html or the more general "pair_style
granular"_pair_granular.html" commands. Currently the options are
{hooke}, {hooke/history}, or {hertz/history} for the former, and
{granular} with all the possible options of the associated
{pair_coeff} command for the latter. The equation for the force
between the wall and particles touching it is the same as the
corresponding equation on the "pair_style gran/*"_pair_gran.html and
"pair_style_granular"_pair_granular.html doc pages, but the effective
radius is calculated using the radius of the particle and the radius
of curvature of the wall at the contact point.
Specifically, delta = radius - r = overlap of particle with wall,
m_eff = mass of particle, and RiRj/Ri+Rj is the effective radius, with
@ -141,12 +151,18 @@ particle.
The parameters {Kn}, {Kt}, {gamma_n}, {gamma_t}, {xmu} and {dampflag}
have the same meaning and units as those specified with the
"pair_style granular"_pair_gran.html commands. This means a NULL can
be used for either {Kt} or {gamma_t} as described on that page. If a
"pair_style gran/*"_pair_gran.html commands. This means a NULL can be
used for either {Kt} or {gamma_t} as described on that page. If a
NULL is used for {Kt}, then a default value is used where {Kt} = 2/7
{Kn}. If a NULL is used for {gamma_t}, then a default value is used
where {gamma_t} = 1/2 {gamma_n}.
All the model choices for cohesion, tangential friction, rolling
friction and twisting friction supported by the "pair_style
granular"_pair_granular.html through its {pair_coeff} command are also
supported for walls. These are discussed in greater detail on the doc
page for "pair_style granular"_pair_granular.html.
Note that you can choose a different force styles and/or different
values for the 6 wall/particle coefficients than for particle/particle
interactions. E.g. if you wish to model the wall as a different
@ -154,9 +170,9 @@ material.
[Restart, fix_modify, output, run start/stop, minimize info:]
Similar to "fix wall/gran"_fix_wall_gran.html command, this fix
writes the shear friction state of atoms interacting with the wall to
"binary restart files"_restart.html, so that a simulation can continue
Similar to "fix wall/gran"_fix_wall_gran.html command, this fix writes
the shear friction state of atoms interacting with the wall to "binary
restart files"_restart.html, so that a simulation can continue
correctly if granular potentials with shear "history" effects are
being used. This fix also includes info about a moving region in the
restart file. See the "read_restart"_read_restart.html command for
@ -170,14 +186,14 @@ So you must re-define your region and if it is a moving region, define
its motion attributes in a way that is consistent with the simulation
that wrote the restart file. In particular, if you want to change the
region motion attributes (e.g. its velocity), then you should ensure
the position/orientation of the region at the initial restart
timestep is the same as it was on the timestep the restart file was
written. If this is not possible, you may need to ignore info in the
restart file by defining a new fix wall/gran/region command in your
restart script, e.g. with a different fix ID. Or if you want to keep
the shear history info but discard the region motion information, you
can use the same fix ID for fix wall/gran/region, but assign it a
region with a different region ID.
the position/orientation of the region at the initial restart timestep
is the same as it was on the timestep the restart file was written.
If this is not possible, you may need to ignore info in the restart
file by defining a new fix wall/gran/region command in your restart
script, e.g. with a different fix ID. Or if you want to keep the
shear history info but discard the region motion information, you can
use the same fix ID for fix wall/gran/region, but assign it a region
with a different region ID.
None of the "fix_modify"_fix_modify.html options are relevant to this
fix. No global or per-atom quantities are stored by this fix for

1
doc/src/lammps.book
View File

@ -580,6 +580,7 @@ pair_extep.html
pair_gauss.html
pair_gayberne.html
pair_gran.html
pair_granular.html
pair_gromacs.html
pair_gw.html
pair_ilp_graphene_hbn.html

765
doc/src/pair_granular.txt Normal file
View File

@ -0,0 +1,765 @@
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
</script>
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Commands_all.html)
:line
pair_style granular command :h3
[Syntax:]
pair_style granular cutoff :pre
cutoff = global cutoff (optional). See discussion below. :l
[Examples:]
pair_style granular
pair_coeff * * hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 :pre
pair_style granular
pair_coeff * * hertz 1000.0 50.0 tangential mindlin NULL 1.0 0.4 :pre
pair_style granular
pair_coeff * * hertz/material 1e8 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji :pre
pair_style granular
pair_coeff 1 1 jkr 1000.0 50.0 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall :pre
pair_style granular
pair_coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 0.5 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall
pair_coeff 1 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin 800.0 0.5 0.1 roll sds 500.0 200.0 0.1 twisting marshall :pre
[Description:]
The {granular} styles support a variety of options for the normal,
tangential, rolling and twisting forces resulting from contact between
two granular particles. This expands on the options offered by the
"pair gran/*"_pair_gran.html pair styles. The total computed forces
and torques are the sum of various models selected for the normal,
tangential, rolling and twisting modes of motion.
All model choices and parameters are entered in the
"pair_coeff"_pair_coeff.html command, as described below. Unlike
e.g. "pair gran/hooke"_pair_gran.html, coefficient values are not
global, but can be set to different values for different combinations
of particle types, as determined by the "pair_coeff"_pair_coeff.html
command. If the contact model choice is the same for two particle
types, the mixing for the cross-coefficients can be carried out
automatically. This is shown in the second example, where model
choices are the same for type 1 - type 1 as for type 2 - type2
interactions, but coefficients are different. In this case, the
coefficients for type 2 - type interactions can be determined from
mixing rules discussed below. For additional flexibility,
coefficients as well as model forms can vary between particle types,
as shown in the third example: type 1- type 1 interactions are based
on a Hertzian normal contact model and 2-2 interactions are based on a
DMT cohesive model (see below). In that example, 1-1 and 2-2
interactions have different model forms, in which case mixing of
coefficients cannot be determined, so 1-2 interactions must be
explicitly defined via the {pair_coeff 1 2} command, otherwise an
error would result.
:line
The first required keyword for the {pair_coeff} command is the normal
contact model. Currently supported options for normal contact models
and their required arguments are:
{hooke} : \(k_n\), \(\eta_\{n0\}\) (or \(e\))
{hertz} : \(k_n\), \(\eta_\{n0\}\) (or \(e\))
{hertz/material} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\)
{dmt} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\), \(\gamma\)
{jkr} : E, \(\eta_\{n0\}\) (or \(e\)), \(\nu\), \(\gamma\) :ol
Here, \(k_n\) is spring stiffness (with units that depend on model
choice, see below); \(\eta_\{n0\}\) is a damping prefactor (or, in its
place a coefficient of restitution \(e\), depending on the choice of
damping mode, see below); E is Young's modulus in units of
{force}/{length}^2, i.e. {pressure}; \(\nu\) is Poisson's ratio and
\(\gamma\) is a surface energy density, in units of
{energy}/{length}^2.
For the {hooke} model, the normal, elastic component of force acting
on particle {i} due to contact with particle {j} is given by:
\begin\{equation\}
\mathbf\{F\}_\{ne, Hooke\} = k_N \delta_\{ij\} \mathbf\{n\}
\end\{equation\}
Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle
overlap, \(R_i, R_j\) are the particle radii, \(\mathbf\{r\}_\{ij\} =
\mathbf\{r\}_i - \mathbf\{r\}_j\) is the vector separating the two
particle centers (note the i-j ordering so that \(F_\{ne\}\) is
positive for repulsion), and \(\mathbf\{n\} =
\frac\{\mathbf\{r\}_\{ij\}\}\{\|\mathbf\{r\}_\{ij\}\|\}\). Therefore,
for {hooke}, the units of the spring constant \(k_n\) are
{force}/{distance}, or equivalently {mass}/{time^2}.
For the {hertz} model, the normal component of force is given by:
\begin\{equation\}
\mathbf\{F\}_\{ne, Hertz\} = k_N R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\}
\end\{equation\}
Here, \(R_\{eff\} = \frac\{R_i R_j\}\{R_i + R_j\}\) is the effective
radius, denoted for simplicity as {R} from here on. For {hertz}, the
units of the spring constant \(k_n\) are {force}/{length}^2, or
equivalently {pressure}.
For the {hertz/material} model, the force is given by:
\begin\{equation\}
\mathbf\{F\}_\{ne, Hertz/material\} = \frac\{4\}\{3\} E_\{eff\} R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\}
\end\{equation\}
Here, \(E_\{eff\} = E = \left(\frac\{1-\nu_i^2\}\{E_i\} +
\frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\}\) is the effective Young's
modulus, with \(\nu_i, \nu_j \) the Poisson ratios of the particles of
types {i} and {j}. Note that if the elastic modulus and the shear
modulus of the two particles are the same, the {hertz/material} model
is equivalent to the {hertz} model with \(k_N = 4/3 E_\{eff\}\)
The {dmt} model corresponds to the
"(Derjaguin-Muller-Toporov)"_#DMT1975 cohesive model, where the force
is simply Hertz with an additional attractive cohesion term:
\begin\{equation\}
\mathbf\{F\}_\{ne, dmt\} = \left(\frac\{4\}\{3\} E R^\{1/2\}\delta_\{ij\}^\{3/2\} - 4\pi\gamma R\right)\mathbf\{n\}
\end\{equation\}
The {jkr} model is the "(Johnson-Kendall-Roberts)"_#JKR1971 model,
where the force is computed as:
\begin\{equation\}
\label\{eq:force_jkr\}
\mathbf\{F\}_\{ne, jkr\} = \left(\frac\{4Ea^3\}\{3R\} - 2\pi a^2\sqrt\{\frac\{4\gamma E\}\{\pi a\}\}\right)\mathbf\{n\}
\end\{equation\}
Here, {a} is the radius of the contact zone, related to the overlap
\(\delta\) according to:
\begin\{equation\}
\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}
\end\{equation\}
LAMMPS internally inverts the equation above to solve for {a} in terms
of \(\delta\), then solves for the force in the previous
equation. Additionally, note that the JKR model allows for a tensile
force beyond contact (i.e. for \(\delta < 0\)), up to a maximum of
\(3\pi\gamma R\) (also known as the 'pull-off' force). Note that this
is a hysteretic effect, where particles that are not contacting
initially will not experience force until they come into contact
\(\delta \geq 0\); as they move apart and (\(\delta < 0\)), they
experience a tensile force up to \(3\pi\gamma R\), at which point they
lose contact.
:line
In addition, the normal force is augmented by a damping term of the
following general form:
\begin\{equation\}
\mathbf\{F\}_\{n,damp\} = -\eta_n \mathbf\{v\}_\{n,rel\}
\end\{equation\}
Here, \(\mathbf\{v\}_\{n,rel\} = (\mathbf\{v\}_j - \mathbf\{v\}_i)
\cdot \mathbf\{n\}\) is the component of relative velocity along
\(\mathbf\{n\}\).
The optional {damping} keyword to the {pair_coeff} command followed by
a keyword determines the model form of the damping factor \(\eta_n\),
and the interpretation of the \(\eta_\{n0\}\) or \(e\) coefficients
specified as part of the normal contact model settings. The {damping}
keyword and corresponding model form selection may be appended
anywhere in the {pair coeff} command. Note that the choice of damping
model affects both the normal and tangential damping (and depending on
other settings, potentially also the twisting damping). The options
for the damping model currently supported are:
{velocity}
{viscoelastic}
{tsuji} :ol
If the {damping} keyword is not specified, the {viscoelastic} model is
used by default.
For {damping velocity}, the normal damping is simply equal to the
user-specified damping coefficient in the {normal} model:
\begin\{equation\}
\eta_n = \eta_\{n0\}\
\end\{equation\}
Here, \(\gamma_n\) is the damping coefficient specified for the normal
contact model, in units of {mass}/{time},
The {damping viscoelastic} model is based on the viscoelastic
treatment of "(Brilliantov et al)"_#Brill1996, where the normal
damping is given by:
\begin\{equation\}
\eta_n = \eta_\{n0\}\ a m_\{eff\}
\end\{equation\}
Here, \(m_\{eff\} = m_i m_j/(m_i + m_j)\) is the effective mass, {a}
is the contact radius, given by \(a =\sqrt\{R\delta\}\) for all models
except {jkr}, for which it is given implicitly according to \(delta =
a^2/R - 2\sqrt\{\pi \gamma a/E\}\). In this case, \eta_\{n0\}\ is in
units of 1/({time}*{distance}).
The {tsuji} model is based on the work of "(Tsuji et
al)"_#Tsuji1992. Here, the damping coefficient specified as part of
the normal model is interpreted as a restitution coefficient
\(e\). The damping constant \(\eta_n\) is given by:
\begin\{equation\}
\eta_n = \alpha (m_\{eff\}k_n)^\{1/2\}
\end\{equation\}
For normal contact models based on material parameters, \(k_n =
4/3Ea\). The parameter \(\alpha\) is related to the restitution
coefficient {e} according to:
\begin\{equation\}
\alpha = 1.2728-4.2783e+11.087e^2-22.348e^3+27.467e^4-18.022e^5+4.8218e^6
\end\{equation\}
The dimensionless coefficient of restitution \(e\) specified as part
of the normal contact model parameters should be between 0 and 1, but
no error check is performed on this.
The total normal force is computed as the sum of the elastic and
damping components:
\begin\{equation\}
\mathbf\{F\}_n = \mathbf\{F\}_\{ne\} + \mathbf\{F\}_\{n,damp\}
\end\{equation\}
:line
The {pair_coeff} command also requires specification of the tangential
contact model. The required keyword {tangential} is expected, followed
by the model choice and associated parameters. Currently supported
tangential model choices and their expected parameters are as follows:
{linear_nohistory} : \(x_\{\gamma,t\}\), \(\mu_s\)
{linear_history} : \(k_t\), \(x_\{\gamma,t\}\), \(\mu_s\)
{mindlin} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\)
{mindlin_rescale} : \(k_t\) or NULL, \(x_\{\gamma,t\}\), \(\mu_s\) :ol
Here, \(x_\{\gamma,t\}\) is a dimensionless multiplier for the normal
damping \(\eta_n\) that determines the magnitude of the tangential
damping, \(\mu_t\) is the tangential (or sliding) friction
coefficient, and \(k_t\) is the tangential stiffness coefficient.
For {tangential linear_nohistory}, a simple velocity-dependent Coulomb
friction criterion is used, which mimics the behavior of the {pair
gran/hooke} style. The tangential force (\mathbf\{F\}_t\) is given by:
\begin\{equation\}
\mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|\mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}
\end\{equation\}
The tangential damping force \(\mathbf\{F\}_\mathrm\{t,damp\}\) is given by:
\begin\{equation\}
\mathbf\{F\}_\mathrm\{t,damp\} = -\eta_t \mathbf\{v\}_\{t,rel\}
\end\{equation\}
The tangential damping prefactor \(\eta_t\) is calculated by scaling
the normal damping \(\eta_n\) (see above):
\begin\{equation\}
\eta_t = -x_\{\gamma,t\} \eta_n
\end\{equation\}
The normal damping prefactor \(\eta_n\) is determined by the choice of
the {damping} keyword, as discussed above. Thus, the {damping}
keyword also affects the tangential damping. The parameter
\(x_\{\gamma,t\}\) is a scaling coefficient. Several works in the
literature use \(x_\{\gamma,t\} = 1\) ("Marshall"_#Marshall2009,
"Tsuji et al"_#Tsuji1992, "Silbert et al"_#Silbert2001). The relative
tangential velocity at the point of contact is given by
\(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i +
R_j\Omega_j) \times \mathbf\{n\}\), where \(\mathbf\{v\}_\{t\} =
\mathbf\{v\}_r - \mathbf\{v\}_r\cdot\mathbf\{n\}\), \(\mathbf\{v\}_r =
\mathbf\{v\}_j - \mathbf\{v\}_i\). The direction of the applied force
is \(\mathbf\{t\} =
\mathbf\{v_\{t,rel\}\}/\|\mathbf\{v_\{t,rel\}\}\|\).
The normal force value \(F_\{n0\}\) used to compute the critical force
depends on the form of the contact model. For non-cohesive models
({hertz}, {hertz/material}, {hooke}), it is given by the magnitude of
the normal force:
\begin\{equation\}
F_\{n0\} = \|\mathbf\{F\}_n\|
\end\{equation\}
For cohesive models such as {jkr} and {dmt}, the critical force is
adjusted so that the critical tangential force approaches \(\mu_t
F_\{pulloff\}\), see "Marshall"_#Marshall2009, equation 43, and
"Thornton"_#Thornton1991. For both models, \(F_\{n0\}\) takes the
form:
\begin\{equation\}
F_\{n0\} = \|\mathbf\{F\}_ne + 2 F_\{pulloff\}\|
\end\{equation\}
Where \(F_\{pulloff\} = 3\pi \gamma R \) for {jkr}, and
\(F_\{pulloff\} = 4\pi \gamma R \) for {dmt}.
The remaining tangential options all use accumulated tangential
displacement (i.e. contact history). This is discussed below in the
context of the {linear_history} option, but the same treatment of the
accumulated displacement applies to the other options as well.
For {tangential linear_history}, the tangential force is given by:
\begin\{equation\}
\mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|-k_t\mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}
\end\{equation\}
Here, \(\mathbf\{\xi\}\) is the tangential displacement accumulated
during the entire duration of the contact:
\begin\{equation\}
\mathbf\{\xi\} = \int_\{t0\}^t \mathbf\{v\}_\{t,rel\}(\tau) \mathrm\{d\}\tau
\end\{equation\}
This accumulated tangential displacement must be adjusted to account
for changes in the frame of reference of the contacting pair of
particles during contact. This occurs due to the overall motion of the
contacting particles in a rigid-body-like fashion during the duration
of the contact. There are two modes of motion that are relevant: the
'tumbling' rotation of the contacting pair, which changes the
orientation of the plane in which tangential displacement occurs; and
'spinning' rotation of the contacting pair about the vector connecting
their centers of mass (\(\mathbf\{n\}\)). Corrections due to the
former mode of motion are made by rotating the accumulated
displacement into the plane that is tangential to the contact vector
at each step, or equivalently removing any component of the tangential
displacement that lies along \(\mathbf\{n\}\), and rescaling to
preserve the magnitude. This follows the discussion in
"Luding"_#Luding2008, see equation 17 and relevant discussion in that
work:
\begin\{equation\}
\mathbf\{\xi\} = \left(\mathbf\{\xi'\} - (\mathbf\{n\} \cdot \mathbf\{\xi'\})\mathbf\{n\}\right) \frac\{\|\mathbf\{\xi'\}\|\}\{\|\mathbf\{\xi'\}\| - \mathbf\{n\}\cdot\mathbf\{\xi'\}\}
\label\{eq:rotate_displacements\}
\end\{equation\}
Here, \(\mathbf\{\xi'\}\) is the accumulated displacement prior to the
current time step and \(\mathbf\{\xi\}\) is the corrected
displacement. Corrections to the displacement due to the second mode
of motion described above (rotations about \(\mathbf\{n\}\)) are not
currently implemented, but are expected to be minor for most
simulations.
Furthermore, when the tangential force exceeds the critical force, the
tangential displacement is re-scaled to match the value for the
critical force (see "Luding"_#Luding2008, equation 20 and related
discussion):
\begin\{equation\}
\mathbf\{\xi\} = -\frac\{1\}\{k_t\}\left(\mu_t F_\{n0\}\mathbf\{t\} + \mathbf\{F\}_\{t,damp\}\right)
\end\{equation\}
The tangential force is added to the total normal force (elastic plus
damping) to produce the total force on the particle. The tangential
force also acts at the contact point (defined as the center of the
overlap region) to induce a torque on each particle according to:
\begin\{equation\}
\mathbf\{\tau\}_i = -(R_i - 0.5 \delta) \mathbf\{n\} \times \mathbf\{F\}_t
\end\{equation\}
\begin\{equation\}
\mathbf\{\tau\}_j = -(R_j - 0.5 \delta) \mathbf\{n\} \times \mathbf\{F\}_t
\end\{equation\}
For {tangential mindlin}, the "Mindlin"_#Mindlin1949 no-slip solution is used, which differs from the {linear_history}
option by an additional factor of {a}, the radius of the contact region. The tangential force is given by:
\begin\{equation\}
\mathbf\{F\}_t = -min(\mu_t F_\{n0\}, \|-k_t a \mathbf\{\xi\} + \mathbf\{F\}_\mathrm\{t,damp\}\|) \mathbf\{t\}
\end\{equation\}
Here, {a} is the radius of the contact region, given by \(a = \delta
R\) for all normal contact models, except for {jkr}, where it is given
implicitly by \(\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\), see
discussion above. To match the Mindlin solution, one should set \(k_t
= 8G\), where \(G\) is the shear modulus, related to Young's modulus
\(E\) by \(G = E/(2(1+\nu))\), where \(\nu\) is Poisson's ratio. This
can also be achieved by specifying {NULL} for \(k_t\), in which case a
normal contact model that specifies material parameters \(E\) and
\(\nu\) is required (e.g. {hertz/material}, {dmt} or {jkr}). In this
case, mixing of the shear modulus for different particle types {i} and
{j} is done according to:
\begin\{equation\}
1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j
\end\{equation\}
The {mindlin_rescale} option uses the same form as {mindlin}, but the
magnitude of the tangential displacement is re-scaled as the contact
unloads, i.e. if \(a < a_\{t_\{n-1\}\}\):
\begin\{equation\}
\mathbf\{\xi\} = \mathbf\{\xi_\{t_\{n-1\}\}\} \frac\{a\}\{a_\{t_\{n-1\}\}\}
\end\{equation\}
Here, \(t_\{n-1\}\) indicates the value at the previous time
step. This rescaling accounts for the fact that a decrease in the
contact area upon unloading leads to the contact being unable to
support the previous tangential loading, and spurious energy is
created without the rescaling above ("Walton"_#WaltonPC ). See also
discussion in "Thornton et al, 2013"_#Thornton2013 , particularly
equation 18(b) of that work and associated discussion.
:line
The optional {rolling} keyword enables rolling friction, which resists
pure rolling motion of particles. The options currently supported are:
{none}
{sds} : \(k_\{roll\}\), \(\gamma_\{roll\}\), \(\mu_\{roll\}\) :ol
If the {rolling} keyword is not specified, the model defaults to {none}.
For {rolling sds}, rolling friction is computed via a
spring-dashpot-slider, using a 'pseudo-force' formulation, as detailed
by "Luding"_#Luding2008. Unlike the formulation in
"Marshall"_#Marshall2009, this allows for the required adjustment of
rolling displacement due to changes in the frame of reference of the
contacting pair. The rolling pseudo-force is computed analogously to
the tangential force:
\begin\{equation\}
\mathbf\{F\}_\{roll,0\} = k_\{roll\} \mathbf\{\xi\}_\{roll\} - \gamma_\{roll\} \mathbf\{v\}_\{roll\}
\end\{equation\}
Here, \(\mathbf\{v\}_\{roll\} = -R(\mathbf\{\Omega\}_i -
\mathbf\{\Omega\}_j) \times \mathbf\{n\}\) is the relative rolling
velocity, as given in "Wang et al"_#Wang2015 and
"Luding"_#Luding2008. This differs from the expressions given by "Kuhn
and Bagi"_#Kuhn2004 and used in "Marshall"_#Marshall2009; see "Wang et
al"_#Wang2015 for details. The rolling displacement is given by:
\begin\{equation\}
\mathbf\{\xi\}_\{roll\} = \int_\{t_0\}^t \mathbf\{v\}_\{roll\} (\tau) \mathrm\{d\} \tau
\end\{equation\}
A Coulomb friction criterion truncates the rolling pseudo-force if it
exceeds a critical value:
\begin\{equation\}
\mathbf\{F\}_\{roll\} = min(\mu_\{roll\} F_\{n,0\}, \|\mathbf\{F\}_\{roll,0\}\|)\mathbf\{k\}
\end\{equation\}
Here, \(\mathbf\{k\} =
\mathbf\{v\}_\{roll\}/\|\mathbf\{v\}_\{roll\}\|\) is the direction of
the pseudo-force. As with tangential displacement, the rolling
displacement is rescaled when the critical force is exceeded, so that
the spring length corresponds the critical force. Additionally, the
displacement is adjusted to account for rotations of the frame of
reference of the two contacting particles in a manner analogous to the
tangential displacement.
The rolling pseudo-force does not contribute to the total force on
either particle (hence 'pseudo'), but acts only to induce an equal and
opposite torque on each particle, according to:
\begin\{equation\}
\tau_\{roll,i\} = R_\{eff\} \mathbf\{n\} \times \mathbf\{F\}_\{roll\}
\end\{equation\}
\begin\{equation\}
\tau_\{roll,j\} = -\tau_\{roll,i\}
\end\{equation\}
:line
The optional {twisting} keyword enables twisting friction, which
resists rotation of two contacting particles about the vector
\(\mathbf\{n\}\) that connects their centers. The options currently
supported are:
{none}
{sds} : \(k_\{twist\}\), \(\gamma_\{twist\}\), \(\mu_\{twist\}\)
{marshall} :ol
If the {twisting} keyword is not specified, the model defaults to {none}.
For both {twisting sds} and {twisting marshall}, a history-dependent
spring-dashpot-slider is used to compute the twisting torque. Because
twisting displacement is a scalar, there is no need to adjust for
changes in the frame of reference due to rotations of the particle
pair. The formulation in "Marshall"_#Marshall2009 therefore provides
the most straightforward treatment:
\begin\{equation\}
\tau_\{twist,0\} = -k_\{twist\}\xi_\{twist\} - \gamma_\{twist\}\Omega_\{twist\}
\end\{equation\}
Here \(\xi_\{twist\} = \int_\{t_0\}^t \Omega_\{twist\} (\tau)
\mathrm\{d\}\tau\) is the twisting angular displacement, and
\(\Omega_\{twist\} = (\mathbf\{\Omega\}_i - \mathbf\{\Omega\}_j) \cdot
\mathbf\{n\}\) is the relative twisting angular velocity. The torque
is then truncated according to:
\begin\{equation\}
\tau_\{twist\} = min(\mu_\{twist\} F_\{n,0\}, \tau_\{twist,0\})
\end\{equation\}
Similar to the sliding and rolling displacement, the angular
displacement is rescaled so that it corresponds to the critical value
if the twisting torque exceeds this critical value:
\begin\{equation\}
\xi_\{twist\} = \frac\{1\}\{k_\{twist\}\} (\mu_\{twist\} F_\{n,0\}sgn(\Omega_\{twist\}) - \gamma_\{twist\}\Omega_\{twist\})
\end\{equation\}
For {twisting sds}, the coefficients \(k_\{twist\}, \gamma_\{twist\}\)
and \(\mu_\{twist\}\) are simply the user input parameters that follow
the {twisting sds} keywords in the {pair_coeff} command.
For {twisting_marshall}, the coefficients are expressed in terms of
sliding friction coefficients, as discussed in
"Marshall"_#Marshall2009 (see equations 32 and 33 of that work):
\begin\{equation\}
k_\{twist\} = 0.5k_ta^2
\end\{equation\}
\begin\{equation\}
\eta_\{twist\} = 0.5\eta_ta^2
\end\{equation\}
\begin\{equation\}
\mu_\{twist\} = \frac\{2\}\{3\}a\mu_t
\end\{equation\}
Finally, the twisting torque on each particle is given by:
\begin\{equation\}
\mathbf\{\tau\}_\{twist,i\} = \tau_\{twist\}\mathbf\{n\}
\end\{equation\}
\begin\{equation\}
\mathbf\{\tau\}_\{twist,j\} = -\mathbf\{\tau\}_\{twist,i\}
\end\{equation\}
:line
LAMMPS automatically sets pairwise cutoff values for {pair_style
granular} based on particle radii (and in the case of {jkr} pull-off
distances). In the vast majority of situations, this is adequate.
However, a cutoff value can optionally be appended to the {pair_style
granular} command to specify a global cutoff (i.e. a cutoff for all
atom types). Additionally, the optional {cutoff} keyword can be passed
to the {pair_coeff} command, followed by a cutoff value. This will
set a pairwise cutoff for the atom types in the {pair_coeff} command.
These options may be useful in some rare cases where the automatic
cutoff determination is not sufficient, e.g. if particle diameters
are being modified via the {fix adapt} command. In that case, the
global cutoff specified as part of the {pair_style granular} command
is applied to all atom types, unless it is overridden for a given atom
type combination by the {cutoff} value specified in the {pair coeff}
command. If {cutoff} is only specified in the {pair coeff} command
and no global cutoff is appended to the {pair_style granular} command,
then LAMMPS will use that cutoff for the specified atom type
combination, and automatically set pairwise cutoffs for the remaining
atom types.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed on the "Speed packages"_Speed_packages.html doc
page. The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.
These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively. They are only enabled if
LAMMPS was built with those packages. See the "Build
package"_Build_package.html doc page for more info.
You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the "-suffix command-line
switch"_Run_options.html when you invoke LAMMPS, or you can use the
"suffix"_suffix.html command in your input script.
See the "Speed packages"_Speed_packages.html doc page for more
instructions on how to use the accelerated styles effectively.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
The "pair_modify"_pair_modify.html mix, shift, table, and tail options
are not relevant for granular pair styles.
Mixing of coefficients is carried out using geometric averaging for
most quantities, e.g. if friction coefficient for type 1-type 1
interactions is set to \(\mu_1\), and friction coefficient for type
2-type 2 interactions is set to \(\mu_2\), the friction coefficient
for type1-type2 interactions is computed as \(\sqrt\{\mu_1\mu_2\}\)
(unless explicitly specified to a different value by a {pair_coeff 1 2
...} command. The exception to this is elastic modulus, only
applicable to {hertz/material}, {dmt} and {jkr} normal contact
models. In that case, the effective elastic modulus is computed as:
\begin\{equation\}
E_\{eff,ij\} = \left(\frac\{1-\nu_i^2\}\{E_i\} + \frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\}
\end\{equation\}
If the {i-j} coefficients \(E_\{ij\}\) and \(\nu_\{ij\}\) are
explicitly specified, the effective modulus is computed as:
\begin\{equation\}
E_\{eff,ij\} = \left(\frac\{1-\nu_\{ij\}^2\}\{E_\{ij\}\} + \frac\{1-\nu_\{ij\}^2\}\{E_\{ij\}\}\right)^\{-1\}
\end\{equation\}
or
\begin\{equation\}
E_\{eff,ij\} = \frac\{E_\{ij\}\}\{2(1-\nu_\{ij\})\}
\end\{equation\}
These pair styles write their information to "binary restart
files"_restart.html, so a pair_style command does not need to be
specified in an input script that reads a restart file.
These pair styles can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. They do not support the
{inner}, {middle}, {outer} keywords.
The single() function of these pair styles returns 0.0 for the energy
of a pairwise interaction, since energy is not conserved in these
dissipative potentials. It also returns only the normal component of
the pairwise interaction force. However, the single() function also
calculates 10 extra pairwise quantities. The first 3 are the
components of the tangential force between particles I and J, acting
on particle I. The 4th is the magnitude of this tangential force.
The next 3 (5-7) are the components of the rolling torque acting on
particle I. The next entry (8) is the magnitude of the rolling torque.
The next entry (9) is the magnitude of the twisting torque acting
about the vector connecting the two particle centers.
The last 3 (10-12) are the components of the vector connecting
the centers of the two particles (x_I - x_J).
These extra quantities can be accessed by the "compute
pair/local"_compute_pair_local.html command, as {p1}, {p2}, ...,
{p12}.
:line
[Restrictions:]
All the granular pair styles are part of the GRANULAR package. It is
only enabled if LAMMPS was built with that package. See the "Build
package"_Build_package.html doc page for more info.
These pair styles require that atoms store torque and angular velocity
(omega) as defined by the "atom_style"_atom_style.html. They also
require a per-particle radius is stored. The {sphere} atom style does
all of this.
This pair style requires you to use the "comm_modify vel
yes"_comm_modify.html command so that velocities are stored by ghost
atoms.
These pair styles will not restart exactly when using the
"read_restart"_read_restart.html command, though they should provide
statistically similar results. This is because the forces they
compute depend on atom velocities. See the
"read_restart"_read_restart.html command for more details.
[Related commands:]
"pair_coeff"_pair_coeff.html
"pair gran/*"_pair_gran.html
[Default:]
For the {pair_coeff} settings: {damping viscoelastic}, {rolling none},
{twisting none}.
[References:]
:link(Brill1996)
[(Brilliantov et al, 1996)] Brilliantov, N. V., Spahn, F., Hertzsch,
J. M., & Poschel, T. (1996). Model for collisions in granular
gases. Physical review E, 53(5), 5382.
:link(Tsuji1992)
[(Tsuji et al, 1992)] Tsuji, Y., Tanaka, T., & Ishida,
T. (1992). Lagrangian numerical simulation of plug flow of
cohesionless particles in a horizontal pipe. Powder technology, 71(3),
239-250.
:link(JKR1971)
[(Johnson et al, 1971)] Johnson, K. L., Kendall, K., & Roberts,
A. D. (1971). Surface energy and the contact of elastic
solids. Proc. R. Soc. Lond. A, 324(1558), 301-313.
:link(DMT1975)
[Derjaguin et al, 1975)] Derjaguin, B. V., Muller, V. M., & Toporov,
Y. P. (1975). Effect of contact deformations on the adhesion of
particles. Journal of Colloid and interface science, 53(2), 314-326.
:link(Luding2008)
[(Luding, 2008)] Luding, S. (2008). Cohesive, frictional powders:
contact models for tension. Granular matter, 10(4), 235.
:link(Marshall2009)
[(Marshall, 2009)] Marshall, J. S. (2009). Discrete-element modeling
of particulate aerosol flows. Journal of Computational Physics,
228(5), 1541-1561.
:link(Silbert2001)
[(Silbert, 2001)] Silbert, L. E., Ertas, D., Grest, G. S., Halsey,
T. C., Levine, D., & Plimpton, S. J. (2001). Granular flow down an
inclined plane: Bagnold scaling and rheology. Physical Review E,
64(5), 051302.
:link(Kuhn2004)
[(Kuhn and Bagi, 2005)] Kuhn, M. R., & Bagi, K. (2004). Contact
rolling and deformation in granular media. International journal of
solids and structures, 41(21), 5793-5820.
:link(Wang2015)
[(Wang et al, 2015)] Wang, Y., Alonso-Marroquin, F., & Guo,
W. W. (2015). Rolling and sliding in 3-D discrete element
models. Particuology, 23, 49-55.
:link(Thornton1991)
[(Thornton, 1991)] Thornton, C. (1991). Interparticle sliding in the
presence of adhesion. J. Phys. D: Appl. Phys. 24 1942
:link(Mindlin1949)
[(Mindlin, 1949)] Mindlin, R. D. (1949). Compliance of elastic bodies
in contact. J. Appl. Mech., ASME 16, 259-268.
:link(Thornton2013)
[(Thornton et al, 2013)] Thornton, C., Cummins, S. J., & Cleary,
P. W. (2013). An investigation of the comparative behaviour of
alternative contact force models during inelastic collisions. Powder
Technology, 233, 30-46.
:link(WaltonPC)
[(Otis R. Walton)] Walton, O.R., Personal Communication

1
doc/src/pairs.txt
View File

@ -42,6 +42,7 @@ Pair Styles :h1
pair_gauss
pair_gayberne
pair_gran
pair_granular
pair_gromacs
pair_gw
pair_hbond_dreiding