Merge branch 'develop' of github.com:lammps/lammps into amd_hip_opt

This commit is contained in:
Stan Gerald Moore
2022-02-21 13:23:51 -07:00
49 changed files with 4817 additions and 466 deletions

View File

@ -5,7 +5,6 @@ mark_as_advanced(OPENCL_LOADER_URL)
mark_as_advanced(OPENCL_LOADER_MD5)
set(INSTALL_LIBOPENCL OFF CACHE BOOL "" FORCE)
set(BUILD_SHARED_LIBS OFF)
include(ExternalCMakeProject)
ExternalCMakeProject(opencl_loader ${OPENCL_LOADER_URL} ${OPENCL_LOADER_MD5} opencl-loader . "")

View File

@ -1,4 +1,4 @@
.TH LAMMPS "1" "7 January 2022" "2022-1-7"
.TH LAMMPS "1" "17 February 2022" "2022-2-17"
.SH NAME
.B LAMMPS
\- Molecular Dynamics Simulator.

View File

@ -124,6 +124,7 @@ OPT.
* :doc:`hbond/dreiding/morse (o) <pair_hbond_dreiding>`
* :doc:`hdnnp <pair_hdnnp>`
* :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>`
* :doc:`ilp/tmd <pair_ilp_tmd>`
* :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>`
* :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>`
* :doc:`lcbop <pair_lcbop>`
@ -241,6 +242,7 @@ OPT.
* :doc:`reaxff (ko) <pair_reaxff>`
* :doc:`rebo (io) <pair_airebo>`
* :doc:`resquared (go) <pair_resquared>`
* :doc:`saip/metal <pair_saip_metal>`
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>`
* :doc:`smd/hertz <pair_smd_hertz>`
* :doc:`smd/tlsph <pair_smd_tlsph>`

View File

@ -11,7 +11,7 @@ of time and requests from the LAMMPS user community.
:maxdepth: 1
Developer_org
Developer_cxx_vs_c_style
Developer_code_design
Developer_parallel
Developer_flow
Developer_write

View File

@ -0,0 +1,433 @@
Code design
-----------
This section explains some of the code design choices in LAMMPS with
the goal of helping developers write new code similar to the existing
code. Please see the section on :doc:`Requirements for contributed
code <Modify_style>` for more specific recommendations and guidelines.
While that section is organized more in the form of a checklist for
code contributors, the focus here is on overall code design strategy,
choices made between possible alternatives, and discussing some
relevant C++ programming language constructs.
Historically, the basic design philosophy of the LAMMPS C++ code was a
"C with classes" style. The motivation was to make it easy to modify
LAMMPS for people without significant training in C++ programming.
Data structures and code constructs were used that resemble the
previous implementation(s) in Fortran. A contributing factor to this
choice also was that at the time, C++ compilers were often not mature
and some of the advanced features contained bugs or did not function
as the standard required. There were also disagreements between
compiler vendors as to how to interpret the C++ standard documents.
However, C++ compilers have now advanced significantly. In 2020 we
decided to to require the C++11 standard as the minimum C++ language
standard for LAMMPS. Since then we have begun to also replace some of
the C-style constructs with equivalent C++ functionality, either from
the C++ standard library or as custom classes or functions, in order
to improve readability of the code and to increase code reuse through
abstraction of commonly used functionality.
.. note::
Please note that as of spring 2022 there is still a sizable chunk
of legacy code in LAMMPS that has not yet been refactored to
reflect these style conventions in full. LAMMPS has a large code
base and many different contributors and there also is a hierarchy
of precedence in which the code is adapted. Highest priority has
been the code in the ``src`` folder, followed by code in packages
in order of their popularity and complexity (simpler code is
adapted sooner), followed by code in the ``lib`` folder. Source
code that is downloaded from external packages or libraries during
compilation is not subject to the conventions discussed here.
Object oriented code
^^^^^^^^^^^^^^^^^^^^
LAMMPS is designed to be an object oriented code. Each simulation is
represented by an instance of the LAMMPS class. When running in
parallel each MPI process creates such an instance. This can be seen
in the ``main.cpp`` file where the core steps of running a LAMMPS
simulation are the following 3 lines of code:
.. code-block:: C++
LAMMPS *lammps = new LAMMPS(argc, argv, lammps_comm);
lammps->input->file();
delete lammps;
The first line creates a LAMMPS class instance and passes the command
line arguments and the global communicator to its constructor. The
second line triggers the LAMMPS instance to process the input (either
from standard input or a provided input file) until the simulation
ends. The third line deletes the LAMMPS instance. The remainder of
the main.cpp file has code for error handling, MPI configuration, and
other special features.
The basic LAMMPS class hierarchy which is created by the LAMMPS class
constructor is shown in :ref:`class-topology`. When input commands
are processed, additional class instances are created, or deleted, or
replaced. Likewise specific member functions of specific classes are
called to trigger actions such creating atoms, computing forces,
computing properties, time-propagating the system, or writing output.
Compositing and Inheritance
===========================
LAMMPS makes extensive use of the object oriented programming (OOP)
principles of *compositing* and *inheritance*. Classes like the
``LAMMPS`` class are a **composite** containing pointers to instances
of other classes like ``Atom``, ``Comm``, ``Force``, ``Neighbor``,
``Modify``, and so on. Each of these classes implement certain
functionality by storing and manipulating data related to the
simulation and providing member functions that trigger certain
actions. Some of those classes like ``Force`` are themselves
composites, containing instances of classes describing different force
interactions. Similarly the ``Modify`` class contains a list of
``Fix`` and ``Compute`` classes. If the input commands that
correspond to these classes include the word *style*, then LAMMPS
stores only a single instance of that class. E.g. *atom_style*,
*comm_style*, *pair_style*, *bond_style*. It the input command does
not include the word *style*, there can be many instances of that
class defined. E.g. *region*, *fix*, *compute*, *dump*.
**Inheritance** enables creation of *derived* classes that can share
common functionality in their base class while providing a consistent
interface. The derived classes replace (dummy or pure) functions in
the base class. The higher level classes can then call those methods
of the instantiated classes without having to know which specific
derived class variant was instantiated. In LAMMPS these derived
classes are often referred to as "styles", e.g. pair styles, fix
styles, atom styles and so on.
This is the origin of the flexibility of LAMMPS. For example pair
styles implement a variety of different non-bonded interatomic
potentials functions. All details for the implementation of a
potential are stored and executed in a single class.
As mentioned above, there can be multiple instances of classes derived
from the ``Fix`` or ``Compute`` base classes. They represent a
different facet of LAMMPS flexibility as they provide methods which
can be called at different points in time within a timestep, as
explained in `Developer_flow`. This allows the input script to tailor
how a specific simulation is run, what diagnostic computations are
performed, and how the output of those computations is further
processed or output.
Additional code sharing is possible by creating derived classes from the
derived classes (e.g., to implement an accelerated version of a pair
style) where only a subset of the derived class methods are replaced
with accelerated versions.
Polymorphism
============
Polymorphism and dynamic dispatch are another OOP feature that play an
important role in how LAMMPS selects what code to execute. In a
nutshell, this is a mechanism where the decision of which member
function to call from a class is determined at runtime and not when
the code is compiled. To enable it, the function has to be declared
as ``virtual`` and all corresponding functions in derived classes
should use the ``override`` property. Below is a brief example.
.. code-block:: c++
class Base {
public:
virtual ~Base() = default;
void call();
void normal();
virtual void poly();
};
void Base::call() {
normal();
poly();
}
class Derived : public Base {
public:
~Derived() override = default;
void normal();
void poly() override;
};
// [....]
Base *base1 = new Base();
Base *base2 = new Derived();
base1->call();
base2->call();
The difference in behavior of the ``normal()`` and the ``poly()`` member
functions is which of the two member functions is called when executing
`base1->call()` versus `base2->call()`. Without polymorphism, a
function within the base class can only call member functions within the
same scope, that is ``Base::call()`` will always call
``Base::normal()``. But for the `base2->call()` case the call of the
virtual member function will be dispatched to ``Derived::poly()``
instead. This mechanism means that functions are called within the
scope of the class type that was used to *create* the class instance are
invoked; even if they are assigned to a pointer using the type of a base
class. This is the desired behavior and this way LAMMPS can even use
styles that are loaded at runtime from a shared object file with the
:doc:`plugin command <plugin>`.
A special case of virtual functions are so-called pure functions. These
are virtual functions that are initialized to 0 in the class declaration
(see example below).
.. code-block:: c++
class Base {
public:
virtual void pure() = 0;
};
This has the effect that an instance of the base class cannot be
created and that derived classes **must** implement these functions.
Many of the functions listed with the various class styles in the
section :doc:`Modify` are pure functions. The motivation for this is
to define the interface or API of the functions but defer their
implementation to the derived classes.
However, there are downsides to this. For example, calls to virtual
functions from within a constructor, will not be in the scope of the
derived class and thus it is good practice to either avoid calling them
or to provide an explicit scope such as ``Base::poly()`` or
``Derived::poly()``. Furthermore, any destructors in classes containing
virtual functions should be declared virtual too, so they will be
processed in the expected order before types are removed from dynamic
dispatch.
.. admonition:: Important Notes
In order to be able to detect incompatibilities at compile time and
to avoid unexpected behavior, it is crucial that all member functions
that are intended to replace a virtual or pure function use the
``override`` property keyword. For the same reason, the use of
overloads or default arguments for virtual functions should be
avoided as they lead to confusion over which function is supposed to
override which and which arguments need to be declared.
Style Factories
===============
In order to create class instances for different styles, LAMMPS often
uses a programming pattern called `Factory`. Those are functions that
create an instance of a specific derived class, say ``PairLJCut`` and
return a pointer to the type of the common base class of that style,
``Pair`` in this case. To associate the factory function with the
style keyword, an ``std::map`` class is used with function pointers
indexed by their keyword (for example "lj/cut" for ``PairLJCut`` and
"morse" for ``PairMorse``). A couple of typedefs help keep the code
readable and a template function is used to implement the actual
factory functions for the individual classes. Below is an example
of such a factory function from the ``Force`` class as declared in
``force.h`` and implemented in ``force.cpp``. The file ``style_pair.h``
is generated during compilation and includes all main header files
(i.e. those starting with ``pair_``) of pair styles and then the
macro ``PairStyle()`` will associate the style name "lj/cut"
with a factory function creating an instance of the ``PairLJCut``
class.
.. code-block:: C++
// from force.h
typedef Pair *(*PairCreator)(LAMMPS *);
typedef std::map<std::string, PairCreator> PairCreatorMap;
PairCreatorMap *pair_map;
// from force.cpp
template <typename S, typename T> static S *style_creator(LAMMPS *lmp)
{
return new T(lmp);
}
// [...]
pair_map = new PairCreatorMap();
#define PAIR_CLASS
#define PairStyle(key, Class) (*pair_map)[#key] = &style_creator<Pair, Class>;
#include "style_pair.h"
#undef PairStyle
#undef PAIR_CLASS
// from pair_lj_cut.h
#ifdef PAIR_CLASS
PairStyle(lj/cut,PairLJCut);
#else
// [...]
Similar code constructs are present in other files like ``modify.cpp`` and
``modify.h`` or ``neighbor.cpp`` and ``neighbor.h``. Those contain
similar macros and include ``style_*.h`` files for creating class instances
of styles they manage.
I/O and output formatting
^^^^^^^^^^^^^^^^^^^^^^^^^
C-style stdio versus C++ style iostreams
========================================
LAMMPS uses the "stdio" library of the standard C library for reading
from and writing to files and console instead of C++ "iostreams".
This is mainly motivated by better performance, better control over
formatting, and less effort to achieve specific formatting.
Since mixing "stdio" and "iostreams" can lead to unexpected
behavior. use of the latter is strongly discouraged. Also output to
the screen should not use the predefined ``stdout`` FILE pointer, but
rather the ``screen`` and ``logfile`` FILE pointers managed by the
LAMMPS class. Furthermore, output should generally only be done by
MPI rank 0 (``comm->me == 0``). Output that is sent to both
``screen`` and ``logfile`` should use the :cpp:func:`utils::logmesg()
convenience function <LAMMPS_NS::utils::logmesg>`.
We also discourage the use of stringstreams because the bundled {fmt}
library and the customized tokenizer classes can provide the same
functionality in a cleaner way with better performance. This also
helps maintain a consistent programming syntax with code from many
different contributors.
Formatting with the {fmt} library
===================================
The LAMMPS source code includes a copy of the `{fmt} library
<https://fmt.dev>`_ which is preferred over formatting with the
"printf()" family of functions. The primary reason is that it allows
a typesafe default format for any type of supported data. This is
particularly useful for formatting integers of a given size (32-bit or
64-bit) which may require different format strings depending on
compile time settings or compilers/operating systems. Furthermore,
{fmt} gives better performance, has more functionality, a familiar
formatting syntax that has similarities to ``format()`` in Python, and
provides a facility that can be used to integrate format strings and a
variable number of arguments into custom functions in a much simpler
way than the varargs mechanism of the C library. Finally, {fmt} has
been included into the C++20 language standard, so changes to adopt it
are future-proof.
Formatted strings are frequently created by calling the
``fmt::format()`` function which will return a string as a
``std::string`` class instance. In contrast to the ``%`` placeholder
in ``printf()``, the {fmt} library uses ``{}`` to embed format
descriptors. In the simplest case, no additional characters are
needed as {fmt} will choose the default format based on the data type
of the argument. Otherwise the ``fmt::print()`` function may be
used instead of ``printf()`` or ``fprintf()``. In addition, several
LAMMPS output functions, that originally accepted a single string as
argument have been overloaded to accept a format string with optional
arguments as well (e.g., ``Error::all()``, ``Error::one()``,
``utils::logmesg()``).
Summary of the {fmt} format syntax
==================================
The syntax of the format string is "{[<argument id>][:<format spec>]}",
where either the argument id or the format spec (separated by a colon
':') is optional. The argument id is usually a number starting from 0
that is the index to the arguments following the format string. By
default these are assigned in order (i.e. 0, 1, 2, 3, 4 etc.). The most
common case for using argument id would be to use the same argument in
multiple places in the format string without having to provide it as an
argument multiple times. In LAMMPS the argument id is rarely used.
More common is the use of a format specifier, which starts with a colon.
This may optionally be followed by a fill character (default is ' '). If
provided, the fill character **must** be followed by an alignment
character ('<', '^', '>' for left, centered, or right alignment
(default)). The alignment character may be used without a fill
character. The next important format parameter would be the minimum
width, which may be followed by a dot '.' and a precision for floating
point numbers. The final character in the format string would be an
indicator for the "presentation", i.e. 'd' for decimal presentation of
integers, 'x' for hexadecimal, 'o' for octal, 'c' for character etc.
This mostly follows the "printf()" scheme but without requiring an
additional length parameter to distinguish between different integer
widths. The {fmt} library will detect those and adapt the formatting
accordingly. For floating point numbers there are correspondingly, 'g'
for generic presentation, 'e' for exponential presentation, and 'f' for
fixed point presentation.
Thus "{:8}" would represent *any* type argument using at least 8
characters; "{:<8}" would do this as left aligned, "{:^8}" as centered,
"{:>8}" as right aligned. If a specific presentation is selected, the
argument type must be compatible or else the {fmt} formatting code will
throw an exception. Some format string examples are given below:
.. code-block:: C
auto mesg = fmt::format(" CPU time: {:4d}:{:02d}:{:02d}\n", cpuh, cpum, cpus);
mesg = fmt::format("{:<8s}| {:<10.5g} | {:<10.5g} | {:<10.5g} |{:6.1f} |{:6.2f}\n",
label, time_min, time, time_max, time_sq, tmp);
utils::logmesg(lmp,"{:>6} = max # of 1-2 neighbors\n",maxall);
utils::logmesg(lmp,"Lattice spacing in x,y,z = {:.8} {:.8} {:.8}\n",
xlattice,ylattice,zlattice);
which will create the following output lines:
.. parsed-literal::
CPU time: 0:02:16
Pair | 2.0133 | 2.0133 | 2.0133 | 0.0 | 84.21
4 = max # of 1-2 neighbors
Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962
Finally, a special feature of the {fmt} library is that format
parameters like the width or the precision may be also provided as
arguments. In that case a nested format is used where a pair of curly
braces (with an optional argument id) "{}" are used instead of the
value, for example "{:{}d}" will consume two integer arguments, the
first will be the value shown and the second the minimum width.
For more details and examples, please consult the `{fmt} syntax
documentation <https://fmt.dev/latest/syntax.html>`_ website.
Memory management
^^^^^^^^^^^^^^^^^
Dynamical allocation of small data and objects can be done with the
the C++ commands "new" and "delete/delete[]. Large data should use
the member functions of the ``Memory`` class, most commonly,
``Memory::create()``, ``Memory::grow()``, and ``Memory::destroy()``,
which provide variants for vectors, 2d arrays, 3d arrays, etc.
These can also be used for small data.
The use of ``malloc()``, ``calloc()``, ``realloc()`` and ``free()``
directly is strongly discouraged. To simplify adapting legacy code
into the LAMMPS code base the member functions ``Memory::smalloc()``,
``Memory::srealloc()``, and ``Memory::sfree()`` are available, which
perform additional error checks for safety.
Use of these custom memory allocation functions is motivated by the
following considerations:
- memory allocation failures on *any* MPI rank during a parallel run
will trigger an immediate abort of the entire parallel calculation
instead of stalling it
- a failing "new" will trigger an exception which is also captured by
LAMMPS and triggers a global abort
- allocation of multi-dimensional arrays will be done in a C compatible
fashion but so that the storage of the actual data is stored in one
large contiguous block. Thus when MPI communication is needed,
the data can be communicated directly (similar to Fortran arrays).
- the "destroy()" and "sfree()" functions may safely be called on NULL
pointers
- the "destroy()" functions will nullify the pointer variables making
"use after free" errors easy to detect
- it is possible to use a larger than default memory alignment (not on
all operating systems, since the allocated storage pointers must be
compatible with ``free()`` for technical reasons)
In the practical implementation of code this means that any pointer
variables that are class members should be initialized to a
``nullptr`` value in their respective constructors. That way it is
safe to call ``Memory::destroy()`` or ``delete[]`` on them before
*any* allocation outside the constructor. This helps prevent memory
leaks.

View File

@ -1,384 +0,0 @@
Code design
-----------
This section discusses some of the code design choices in LAMMPS and
overall strategy in order to assist developers to write new code that
will fit well with the remaining code. Please see the section on
:doc:`Requirements for contributed code <Modify_style>` for more
specific recommendations and guidelines. While that section is
organized more in the form of a checklist for code contributors, the
focus here is on overall code design strategy, choices made between
possible alternatives, and to discuss of some relevant C++ programming
language constructs.
Historically, the basic design philosophy of the LAMMPS C++ code was
that of a "C with classes" style. The was motivated by the desire to
make it easier to modify LAMMPS for people without significant training
in C++ programming and by trying to use data structures and code constructs
that somewhat resemble the previous implementation(s) in Fortran.
A contributing factor for this choice also was that at the time the
implementation of C++ compilers was not always very mature and some of
the advanced features contained bugs or were not functioning exactly
as the standard required; plus there was some disagreement between
compiler vendors about how to interpret the C++ standard documents.
However, C++ compilers have advanced a lot since then and with the
transition to requiring the C++11 standard in 2020 as the minimum C++ language
standard for LAMMPS, the decision was made to also replace some of the
C-style constructs with equivalent C++ functionality, either from the
C++ standard library or as custom classes or function, in order to
improve readability of the code and to increase code reuse through
abstraction of commonly used functionality.
.. note::
Please note that as of spring 2022 there is still a sizable chunk of
legacy code in LAMMPS that has not yet been refactored to reflect these
style conventions in full. LAMMPS has a large code base and many
different contributors and there also is a hierarchy of precedence
in which the code is adapted. Highest priority has the code in the
``src`` folder, followed by code in packages in order of their popularity
and complexity (simpler code is adapted sooner), followed by code
in the ``lib`` folder. Source code that is downloaded during compilation
is not subject to the conventions discussed here.
Object oriented code
^^^^^^^^^^^^^^^^^^^^
LAMMPS is designed to be an object oriented code, that is each
simulation is represented by an instance of the LAMMPS class. When
running in parallel, of course, each MPI process will create such an
instance. This can be seen in the ``main.cpp`` file where the core
steps of running a LAMMPS simulation are the following 3 lines of code:
.. code-block:: C++
LAMMPS *lammps = new LAMMPS(argc, argv, lammps_comm);
lammps->input->file();
delete lammps;
The first line creates a LAMMPS class instance and passes the command
line arguments and the global communicator to its constructor. The
second line tells the LAMMPS instance to process the input (either from
standard input or the provided input file) until the end. And the third
line deletes that instance again. The remainder of the main.cpp file
are for error handling, MPI configuration and other special features.
In the constructor of the LAMMPS class instance the basic LAMMPS class hierarchy
is created as shown in :ref:`class-topology`. While processing the input further
class instances are created, or deleted, or replaced and specific member functions
of specific classes are called to trigger actions like creating atoms, computing
forces, computing properties, propagating the system, or writing output.
Compositing and Inheritance
===========================
LAMMPS makes extensive use of the object oriented programming (OOP)
principles of *compositing* and *inheritance*. Classes like the
``LAMMPS`` class are a **composite** containing pointers to instances of
other classes like ``Atom``, ``Comm``, ``Force``, ``Neighbor``,
``Modify``, and so on. Each of these classes implement certain
functionality by storing and manipulating data related to the simulation
and providing member functions that trigger certain actions. Some of
those classes like ``Force`` are a composite again containing instances
of classes describing the force interactions or ``Modify`` containing
and calling fixes and computes. In most cases (e.g. ``AtomVec``, ``Comm``,
``Pair``, or ``Bond``) there is only one instance of those member classes
allowed, but in a few cases (e.g. ``Region``, ``Fix``, ``Compute``, or
``Dump``) there can be multiple instances and the parent class is
maintaining a list of the pointers of instantiated classes instead
of a single pointer.
Changing behavior or adjusting how LAMMPS handles a simulation is
implemented via **inheritance** where different variants of the
functionality are realized by creating *derived* classes that can share
common functionality in their base class and provide a consistent
interface where the derived classes replace (dummy or pure) functions in
the base class. The higher level classes can then call those methods of
the instantiated classes without having to know which specific derived
class variant was instantiated. In the LAMMPS documentation those
derived classes are usually referred to a "styles", e.g. pair styles,
fix styles, atom styles and so on.
This is the origin of the flexibility of LAMMPS and facilitates for
example to compute forces for very different non-bonded potential
functions by having different pair styles (implemented as different
classes derived from the ``Pair`` class) where the evaluation of the
potential function is confined to the implementation of the individual
classes. Whenever a new :doc:`pair_style` or :doc:`bond_style` or
:doc:`comm_style` or similar command is processed in the LAMMPS input
any existing class instance is deleted and a new instance created in
it place.
Classes derived from ``Fix`` or ``Compute`` represent a different facet
of LAMMPS' flexibility as there can be multiple instances of them an
their member functions will be called at different phases of the time
integration process (as explained in `Developer_flow`). This way
multiple manipulations of the entire or parts of the system can be
programmed (with fix styles) or different computations can be performed
and accessed and further processed or output through a common interface
(with compute styles).
Further code sharing is possible by creating derived classes from the
derived classes (for instance to implement an accelerated version of a
pair style) where then only a subset of the methods are replaced with
the accelerated versions.
Polymorphism
============
Polymorphism and dynamic dispatch are another OOP feature that play an
important part of how LAMMPS selects which code to execute. In a nutshell,
this is a mechanism where the decision of which member function to call
from a class is determined at runtime and not when the code is compiled.
To enable it, the function has to be declared as ``virtual`` and all
corresponding functions in derived classes should be using the ``override``
property. Below is a brief example.
.. code-block:: c++
class Base {
public:
virtual ~Base() = default;
void call();
void normal();
virtual void poly();
};
void Base::call() {
normal();
poly();
}
class Derived : public Base {
public:
~Derived() override = default;
void normal();
void poly() override;
};
// [....]
Base *base1 = new Base();
Base *base2 = new Derived();
base1->call();
base2->call();
The difference in behavior of the ``normal()`` and the ``poly()`` member
functions is in which of the two member functions is called when
executing `base1->call()` and `base2->call()`. Without polymorphism, a
function within the base class will call only member functions within
the same scope, that is ``Base::call()`` will always call
``Base::normal()``. But for the `base2->call()` the call for the
virtual member function will be dispatched to ``Derived::poly()``
instead. This mechanism allows to always call functions within the
scope of the class type that was used to create the class instance, even
if they are assigned to a pointer using the type of a base class. This
is the desired behavior, and thanks to dynamic dispatch, LAMMPS can even
use styles that are loaded at runtime from a shared object file with the
:doc:`plugin command <plugin>`.
A special case of virtual functions are so-called pure functions. These
are virtual functions that are initialized to 0 in the class declaration
(see example below).
.. code-block:: c++
class Base {
public:
virtual void pure() = 0;
};
This has the effect that it will no longer be possible to create an
instance of the base class and that derived classes **must** implement
these functions. Many of the functions listed with the various class
styles in the section :doc:`Modify` are such pure functions. The
motivation for this is to define the interface or API of the functions
but defer the implementation to the derived classes.
However, there are downsides to this. For example, calls to virtual
functions from within a constructor, will not be in the scope of the
derived class and thus it is good practice to either avoid calling them
or to provide an explicit scope like in ``Base::poly()``. Furthermore,
any destructors in classes containing virtual functions should be
declared virtual, too, so they are processed in the expected order
before types are removed from dynamic dispatch.
.. admonition:: Important Notes
In order to be able to detect incompatibilities and to avoid unexpected
behavior already at compile time, it is crucial that all member functions
that are intended to replace a virtual or pure function use the ``override``
property keyword. For the same reason it should be avoided to use overloads
or default arguments for virtual functions as they lead to confusion over
which function is supposed to override which and which arguments need to be
declared.
Style Factories
===============
In order to create class instances of the different styles, LAMMPS often
uses a programming pattern called `Factory`. Those are functions that create
an instance of a specific derived class, say ``PairLJCut`` and return a pointer
to the type of the common base class of that style, ``Pair`` in this case.
To associate the factory function with the style keyword, an ``std::map``
class is used in which function pointers are indexed by their keyword
(for example "lj/cut" for ``PairLJCut`` and "morse" ``PairMorse``).
A couple of typedefs help to keep the code readable and a template function
is used to implement the actual factory functions for the individual classes.
I/O and output formatting
^^^^^^^^^^^^^^^^^^^^^^^^^
C-style stdio versus C++ style iostreams
========================================
LAMMPS chooses to use the "stdio" library of the standard C library for
reading from and writing to files and console instead of C++
"iostreams". This is mainly motivated by the better performance, better
control over formatting, and less effort to achieve specific formatting.
Since mixing "stdio" and "iostreams" can lead to unexpected behavior using
the latter is strongly discouraged. Also output to the screen should not
use the predefined ``stdout`` FILE pointer, but rather the ``screen`` and
``logfile`` FILE pointers managed by the LAMMPS class. Furthermore, output
should only be done by MPI rank 0 (``comm->me == 0``) and output that is
send to both ``screen`` and ``logfile`` should use the
:cpp:func:`utils::logmesg() convenience function <LAMMPS_NS::utils::logmesg>`.
We also discourage the use for stringstreams as the bundled {fmt} library
and the customized tokenizer classes can provide the same functionality
in a cleaner way with better performance. This will also help to retain
a consistent programming style despite the many different contributors.
Formatting with the {fmt} library
===================================
The LAMMPS source code includes a copy of the `{fmt} library
<https://fmt.dev>`_ which is preferred over formatting with the
"printf()" family of functions. The primary reason is that it allows a
typesafe default format for any type of supported data. This is
particularly useful for formatting integers of a given size (32-bit or
64-bit) which may require different format strings depending on compile
time settings or compilers/operating systems. Furthermore, {fmt} gives
better performance, has more functionality, a familiar formatting syntax
that has similarities to ``format()`` in Python, and provides a facility
that can be used to integrate format strings and a variable number of
arguments into custom functions in a much simpler way that the varargs
mechanism of the C library. Finally, {fmt} has been included into the
C++20 language standard, so changes to adopt it are future proof.
Formatted strings are frequently created by calling the
``fmt::format()`` function which will return a string as ``std::string``
class instance. In contrast to the ``%`` placeholder in ``printf()``,
the {fmt} library uses ``{}`` to embed format descriptors. In the
simplest case, no additional characters are needed as {fmt} will choose
the default format based on the data type of the argument. Alternatively
The ``fmt::print()`` function may be used instead of ``printf()`` or
``fprintf()``. In addition, several LAMMPS output functions, that
originally accepted a single string as arguments have been overloaded to
accept a format string with optional arguments as well (e.g.
``Error::all()``, ``Error::one()``, ``utils::logmesg()``).
Summary of the {fmt} format syntax
==================================
The syntax of the format string is "{[<argument id>][:<format spec>]}",
where either the argument id or the format spec (separated by a colon
':') is optional. The argument id is usually a number starting from 0
that is the index to the arguments following the format string. By
default these are assigned in order (i.e. 0, 1, 2, 3, 4 etc.). The most
common case for using argument id would be to use the same argument in
multiple places in the format string without having to provide it as an
argument multiple times. In LAMMPS the argument id is rarely used.
More common is the use of the format specifier, which starts with a
colon. This may optionally be followed by a fill character (default is
' '). If provided, the fill character **must** be followed by an
alignment character ('<', '^', '>' for left, centered, or right
alignment (default)). The alignment character may be used without a fill
character. The next important format parameter would be the minimum
width, which may be followed by a dot '.' and a precision for floating
point numbers. The final character in the format string would be an
indicator for the "presentation", i.e. 'd' for decimal presentation of
integers, 'x' for hexadecimal, 'o' for octal, 'c' for character
etc. This mostly follows the "printf()" scheme but without requiring an
additional length parameter to distinguish between different integer
widths. The {fmt} library will detect those and adapt the formatting
accordingly. For floating point numbers there are correspondingly, 'g'
for generic presentation, 'e' for exponential presentation, and 'f' for
fixed point presentation.
Thus "{:8}" would represent *any* type argument using at least 8
characters; "{:<8}" would do this as left aligned, "{:^8}" as centered,
"{:>8}" as right aligned. If a specific presentation is selected, the
argument type must be compatible or else the {fmt} formatting code will
throw an exception. Some format string examples are given below:
.. code-block:: C
auto mesg = fmt::format(" CPU time: {:4d}:{:02d}:{:02d}\n", cpuh, cpum, cpus);
mesg = fmt::format("{:<8s}| {:<10.5g} | {:<10.5g} | {:<10.5g} |{:6.1f} |{:6.2f}\n",
label, time_min, time, time_max, time_sq, tmp);
utils::logmesg(lmp,"{:>6} = max # of 1-2 neighbors\n",maxall);
utils::logmesg(lmp,"Lattice spacing in x,y,z = {:.8} {:.8} {:.8}\n",
xlattice,ylattice,zlattice);
which will create the following output lines:
.. parsed-literal::
CPU time: 0:02:16
Pair | 2.0133 | 2.0133 | 2.0133 | 0.0 | 84.21
4 = max # of 1-2 neighbors
Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962
A special feature of the {fmt} library is that format parameters like
the width or the precision may be also provided as arguments. In that
case a nested format is used where a pair of curly braces (with an
optional argument id) "{}" are used instead of the value, for example
"{:{}d}" will consume two integer arguments, the first will be the value
shown and the second the minimum width.
For more details and examples, please consult the `{fmt} syntax
documentation <https://fmt.dev/latest/syntax.html>`_ website.
Memory management
^^^^^^^^^^^^^^^^^
Dynamical allocation of data and objects should be done with either the
C++ commands "new" and "delete/delete[]" or using member functions of
the ``Memory`` class, most commonly, ``Memory::create()``,
``Memory::grow()``, and ``Memory::destroy()``. The use of ``malloc()``,
``calloc()``, ``realloc()`` and ``free()`` directly is strongly
discouraged. To simplify adapting legacy code into the LAMMPS code base
the member functions ``Memory::smalloc()``, ``Memory::srealloc()``, and
``Memory::sfree()`` are available.
Using those custom memory allocation functions is motivated by the
following considerations:
- memory allocation failures on *any* MPI rank during a parallel run
will trigger an immediate abort of the entire parallel calculation
instead of stalling it
- a failing "new" will trigger an exception which is also captured by
LAMMPS and triggers a global abort
- allocation of multi-dimensional arrays will be done in a C compatible
fashion but so that the storage of the actual data is stored in one
large consecutive block and thus when MPI communication is needed,
only this storage needs to be communicated (similar to Fortran arrays)
- the "destroy()" and "sfree()" functions may safely be called on NULL
pointers
- the "destroy()" functions will nullify the pointer variables making
"use after free" errors easy to detect
- it is possible to use a larger than default memory alignment (not on
all operating systems, since the allocated storage pointers must be
compatible with ``free()`` for technical reasons)
In the practical implementation of code this means that any pointer variables
that are class members should be initialized to a ``nullptr`` value in their
respective constructors. That way it would be safe to call ``Memory::destroy()``
or ``delete[]`` on them before *any* allocation outside the constructor.
This helps to prevent memory leaks.

View File

@ -11,9 +11,9 @@ Reading and parsing of text and text files
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
It is frequently required for a class in LAMMPS to read in additional
data from a file, most commonly potential parameters from a potential
file for manybody potentials. LAMMPS provides several custom classes
and convenience functions to simplify the process. This offers the
data from a file, e.g. potential parameters from a potential file for
manybody potentials. LAMMPS provides several custom classes and
convenience functions to simplify the process. They offer the
following benefits:
- better code reuse and fewer lines of code needed to implement reading
@ -23,24 +23,26 @@ following benefits:
text to a number or returning a 0 on unrecognized text and thus reading incorrect values
- re-entrant code through avoiding global static variables (as used by ``strtok()``)
- transparent support for translating unsupported UTF-8 characters to their ASCII equivalents
(the text to value conversion functions **only** accept ASCII characters)
(the text-to-value conversion functions **only** accept ASCII characters)
In most cases (e.g. potential files) the same data is needed on all MPI
ranks. Then it is best to do the reading and parsing only on MPI rank
0, and communicate the data later with one or more ``MPI_Bcast()``
calls. For reading generic text and potential parameter files the
custom classes :cpp:class:`TextFileReader <LAMMPS_NS::TextFileReader>`
and :cpp:class:`PotentialFileReader <LAMMPS_NS::PotentialFileReader>`
are available. Those classes allow to read the file as individual lines
for which they can return a tokenizer class (see below) for parsing the
line, or they can return blocks of numbers as a vector directly. The
documentation on `File reader classes <file-reader-classes>`_ contains
an example for a typical case.
In most cases (e.g. potential files) the same data is needed on all
MPI ranks. Then it is best to do the reading and parsing only on MPI
rank 0, and communicate the data later with one or more
``MPI_Bcast()`` calls. For reading generic text and potential
parameter files the custom classes :cpp:class:`TextFileReader
<LAMMPS_NS::TextFileReader>` and :cpp:class:`PotentialFileReader
<LAMMPS_NS::PotentialFileReader>` are available. They allow reading
the file as individual lines for which they can return a tokenizer
class (see below) for parsing the line. Or they can return blocks of
numbers as a vector directly. The documentation on `File reader
classes <file-reader-classes>`_ contains an example for a typical
case.
When reading per-atom data, the data in the file usually needs include
an atom ID so it can be associated with a particular atom. In that case
the data can be read in multi-line chunks and broadcast to all MPI ranks
with :cpp:func:`utils::read_lines_from_file()
When reading per-atom data, the data on each line of the file usually
needs to include an atom ID so it can be associated with a particular
atom. In that case the data can be read in multi-line chunks and
broadcast to all MPI ranks with
:cpp:func:`utils::read_lines_from_file()
<LAMMPS_NS::utils::read_lines_from_file>`. Those chunks are then
split into lines, parsed, and applied only to atoms the MPI rank
"owns".
@ -49,15 +51,16 @@ For splitting a string (incrementally) into words and optionally
converting those to numbers, the :cpp:class:`Tokenizer
<LAMMPS_NS::Tokenizer>` and :cpp:class:`ValueTokenizer
<LAMMPS_NS::ValueTokenizer>` can be used. Those provide a superset of
the functionality of ``strtok()`` from the C-library and the latter also
includes conversion to different types. Any errors while processing the
string in those classes will result in an exception, which can be caught
and the error processed as needed. Unlike the C-library functions
``atoi()``, ``atof()``, ``strtol()``, or ``strtod()`` the conversion
will check if the converted text is a valid integer of floating point
number and will not silently return an unexpected or incorrect value.
For example, ``atoi()`` will return 12 when converting "12.5" while the
ValueTokenizer class will throw an :cpp:class:`InvalidIntegerException
the functionality of ``strtok()`` from the C-library and the latter
also includes conversion to different types. Any errors while
processing the string in those classes will result in an exception,
which can be caught and the error processed as needed. Unlike the
C-library functions ``atoi()``, ``atof()``, ``strtol()``, or
``strtod()`` the conversion will check if the converted text is a
valid integer or floating point number and will not silently return an
unexpected or incorrect value. For example, ``atoi()`` will return 12
when converting "12.5", while the ValueTokenizer class will throw an
:cpp:class:`InvalidIntegerException
<LAMMPS_NS::InvalidIntegerException>` if
:cpp:func:`ValueTokenizer::next_int()
<LAMMPS_NS::ValueTokenizer::next_int>` is called on the same string.

View File

@ -159,6 +159,8 @@ Related commands
:doc:`pair_none <pair_none>`,
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
:doc:`pair_style drip <pair_drip>`,
:doc:`pair_style ilp_tmd <pair_ilp_tmd>`,
:doc:`pair_style saip_metal <pair_saip_metal>`,
:doc:`pair_style pair_kolmogorov_crespi_z <pair_kolmogorov_crespi_z>`,
:doc:`pair_style pair_kolmogorov_crespi_full <pair_kolmogorov_crespi_full>`,
:doc:`pair_style pair_lebedeva_z <pair_lebedeva_z>`,

157
doc/src/pair_ilp_tmd.rst Normal file
View File

@ -0,0 +1,157 @@
.. index:: pair_style ilp/tmd
pair_style ilp/tmd command
===================================
Syntax
""""""
.. code-block:: LAMMPS
pair_style [hybrid/overlay ...] ilp/tmd cutoff tap_flag
* cutoff = global cutoff (distance units)
* tap_flag = 0/1 to turn off/on the taper function
Examples
""""""""
.. code-block:: LAMMPS
pair_style hybrid/overlay ilp/tmd 16.0 1
pair_coeff * * ilp/tmd TMD.ILP Mo S S
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
pair_coeff * * ilp/tmd TMD.ILP Mo S S Mo S S
Description
"""""""""""
The *ilp/tmd* style computes the registry-dependent interlayer
potential (ILP) potential for transition metal dichalcogenides (TMD)
as described in :ref:`(Ouyang7) <Ouyang7>`.
.. math::
E = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\
V_{ij} = & {\rm Tap}(r_{ij})\left \{ e^{-\alpha (r_{ij}/\beta -1)}
\left [ \epsilon + f(\rho_{ij}) + f(\rho_{ji})\right ] -
\frac{1}{1+e^{-d\left [ \left ( r_{ij}/\left (s_R \cdot r^{eff} \right ) \right )-1 \right ]}}
\cdot \frac{C_6}{r^6_{ij}} \right \}\\
\rho_{ij}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_i)^2 \\
\rho_{ji}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_j)^2 \\
f(\rho) = & C e^{ -( \rho / \delta )^2 } \\
{\rm Tap}(r_{ij}) = & 20\left ( \frac{r_{ij}}{R_{cut}} \right )^7 -
70\left ( \frac{r_{ij}}{R_{cut}} \right )^6 +
84\left ( \frac{r_{ij}}{R_{cut}} \right )^5 -
35\left ( \frac{r_{ij}}{R_{cut}} \right )^4 + 1
Where :math:`\mathrm{Tap}(r_{ij})` is the taper function which provides
a continuous cutoff (up to third derivative) for interatomic separations
larger than :math:`r_c` :doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`.
It is important to include all the pairs to build the neighbor list for
calculating the normals.
.. note::
Since each MX2 (M = Mo, W and X = S, Se Te) layer contains two
sub-layers of X atoms and one sub-layer of M atoms, the definition of the
normal vectors used for graphene and h-BN is no longer valid for TMDs.
In :ref:`(Ouyang7) <Ouyang7>`, a new definition is proposed, where for
each atom `i`, its six nearest neighboring atoms belonging to the same
sub-layer are chosen to define the normal vector `{\bf n}_i`.
The parameter file (e.g. TMD.ILP), is intended for use with *metal*
:doc:`units <units>`, with energies in meV. Two additional parameters,
*S*, and *rcut* are included in the parameter file. *S* is designed to
facilitate scaling of energies. *rcut* is designed to build the neighbor
list for calculating the normals for each atom pair.
.. note::
The parameters presented in the parameter file (e.g. TMD.ILP),
are fitted with taper function by setting the cutoff equal to 16.0
Angstrom. Using different cutoff or taper function should be careful.
These parameters provide a good description in both short- and long-range
interaction regimes. This feature is essential for simulations in high pressure
regime (i.e., the interlayer distance is smaller than the equilibrium
distance). The benchmark tests and comparison of these parameters can
be found in :ref:`(Ouyang7) <Ouyang7>`.
This potential must be used in combination with hybrid/overlay.
Other interactions can be set to zero using pair_style *none*\ .
This pair style tallies a breakdown of the total interlayer potential
energy into sub-categories, which can be accessed via the :doc:`compute pair <compute_pair>` command as a vector of values of length 2.
The 2 values correspond to the following sub-categories:
1. *E_vdW* = vdW (attractive) energy
2. *E_Rep* = Repulsive energy
To print these quantities to the log file (with descriptive column
headings) the following commands could be included in an input script:
.. code-block:: LAMMPS
compute 0 all pair ilp/tmd
variable Evdw equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step temp epair v_Erep v_Evdw
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair style does not support the pair_modify mix, shift, table, and
tail options.
This pair style does not write their information to binary restart
files, since it is stored in potential files. Thus, you need to
re-specify the pair_style and pair_coeff commands in an input script
that reads a restart file.
Restrictions
""""""""""""
This pair style is part of the INTERLAYER package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair style requires the newton setting to be *on* for pair
interactions.
The TMD.ILP potential file provided with LAMMPS (see the potentials
directory) are parameterized for *metal* units. You can use this
potential with any LAMMPS units, but you would need to create your
BNCH.ILP potential file with coefficients listed in the appropriate
units, if your simulation does not use *metal* units.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`,
:doc:`pair_none <pair_none>`,
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
:doc:`pair_style drip <pair_drip>`,
:doc:`pair_style saip_metal <pair_saip_metal>`,
:doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`,
:doc:`pair_style pair_kolmogorov_crespi_z <pair_kolmogorov_crespi_z>`,
:doc:`pair_style pair_kolmogorov_crespi_full <pair_kolmogorov_crespi_full>`,
:doc:`pair_style pair_lebedeva_z <pair_lebedeva_z>`,
:doc:`pair_style pair_coul_shield <pair_coul_shield>`.
Default
"""""""
tap_flag = 1
----------
.. _Ouyang7:
**(Ouyang7)** W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).

156
doc/src/pair_saip_metal.rst Normal file
View File

@ -0,0 +1,156 @@
.. index:: pair_style saip/metal
pair_style saip/metal command
===================================
Syntax
""""""
.. code-block:: LAMMPS
pair_style [hybrid/overlay ...] saip/metal cutoff tap_flag
* cutoff = global cutoff (distance units)
* tap_flag = 0/1 to turn off/on the taper function
Examples
""""""""
.. code-block:: LAMMPS
pair_style hybrid/overlay saip/metal 16.0 1
pair_coeff * * saip/metal CHAu.ILP Au C H
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam Au_u3.eam Au NULL NULL
pair_coeff * * rebo CH.rebo NULL C H
pair_coeff * * saip/metal CHAu.ILP Au C H
Description
"""""""""""
The *saip/metal* style computes the registry-dependent interlayer
potential (ILP) potential for hetero-junctions formed with hexagonal
2D materials and metal surfaces, as described in :ref:`(Ouyang6) <Ouyang6>`.
.. math::
E = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\
V_{ij} = & {\rm Tap}(r_{ij})\left \{ e^{-\alpha (r_{ij}/\beta -1)}
\left [ \epsilon + f(\rho_{ij}) + f(\rho_{ji})\right ] -
\frac{1}{1+e^{-d\left [ \left ( r_{ij}/\left (s_R \cdot r^{eff} \right ) \right )-1 \right ]}}
\cdot \frac{C_6}{r^6_{ij}} \right \}\\
\rho_{ij}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_i)^2 \\
\rho_{ji}^2 = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_j)^2 \\
f(\rho) = & C e^{ -( \rho / \delta )^2 } \\
{\rm Tap}(r_{ij}) = & 20\left ( \frac{r_{ij}}{R_{cut}} \right )^7 -
70\left ( \frac{r_{ij}}{R_{cut}} \right )^6 +
84\left ( \frac{r_{ij}}{R_{cut}} \right )^5 -
35\left ( \frac{r_{ij}}{R_{cut}} \right )^4 + 1
Where :math:`\mathrm{Tap}(r_{ij})` is the taper function which provides
a continuous cutoff (up to third derivative) for interatomic separations
larger than :math:`r_c` :doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`.
It is important to include all the pairs to build the neighbor list for
calculating the normals.
.. note::
To account for the isotropic nature of the isolated gold atom
electron cloud, their corresponding normal vectors (`{\bf n}_i`) are
assumed to lie along the interatomic vector `{\bf r}_ij`. Notably, this
assumption is suitable for many bulk material surfaces, for
example, for systems possessing s-type valence orbitals or
metallic surfaces, whose valence electrons are mostly
delocalized, such that their Pauli repulsion with the electrons
of adjacent surfaces are isotropic. Caution should be used in
the case of very small gold contacts, for example, nano-clusters,
where edge effects may become relevant.
The parameter file (e.g. CHAu.ILP), is intended for use with *metal*
:doc:`units <units>`, with energies in meV. Two additional parameters,
*S*, and *rcut* are included in the parameter file. *S* is designed to
facilitate scaling of energies. *rcut* is designed to build the neighbor
list for calculating the normals for each atom pair.
.. note::
The parameters presented in the parameter file (e.g. BNCH.ILP),
are fitted with taper function by setting the cutoff equal to 16.0
Angstrom. Using different cutoff or taper function should be careful.
This potential must be used in combination with hybrid/overlay.
Other interactions can be set to zero using pair_style *none*\ .
This pair style tallies a breakdown of the total interlayer potential
energy into sub-categories, which can be accessed via the :doc:`compute pair <compute_pair>` command as a vector of values of length 2.
The 2 values correspond to the following sub-categories:
1. *E_vdW* = vdW (attractive) energy
2. *E_Rep* = Repulsive energy
To print these quantities to the log file (with descriptive column
headings) the following commands could be included in an input script:
.. code-block:: LAMMPS
compute 0 all pair saip/metal
variable Evdw equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step temp epair v_Erep v_Evdw
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair style does not support the pair_modify mix, shift, table, and
tail options.
This pair style does not write their information to binary restart
files, since it is stored in potential files. Thus, you need to
re-specify the pair_style and pair_coeff commands in an input script
that reads a restart file.
Restrictions
""""""""""""
This pair style is part of the INTERLAYER package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair style requires the newton setting to be *on* for pair
interactions.
The CHAu.ILP potential file provided with LAMMPS (see the potentials
directory) are parameterized for *metal* units. You can use this
potential with any LAMMPS units, but you would need to create your
BNCH.ILP potential file with coefficients listed in the appropriate
units, if your simulation does not use *metal* units.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`,
:doc:`pair_none <pair_none>`,
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
:doc:`pair_style drip <pair_drip>`,
:doc:`pair_style ilp_tmd <pair_ilp_tmd>`,
:doc:`pair_style ilp_graphene_hbn <pair_ilp_graphene_hbn>`,
:doc:`pair_style pair_kolmogorov_crespi_z <pair_kolmogorov_crespi_z>`,
:doc:`pair_style pair_kolmogorov_crespi_full <pair_kolmogorov_crespi_full>`,
:doc:`pair_style pair_lebedeva_z <pair_lebedeva_z>`,
:doc:`pair_style pair_coul_shield <pair_coul_shield>`.
Default
"""""""
tap_flag = 1
----------
.. _Ouyang6:
**(Ouyang6)** W. Ouyang, O. Hod, and R. Guerra, J. Chem. Theory Comput. 17, 7215 (2021).

View File

@ -188,6 +188,7 @@ accelerated styles exist.
* :doc:`hbond/dreiding/morse <pair_hbond_dreiding>` - DREIDING hydrogen bonding Morse potential
* :doc:`hdnnp <pair_hdnnp>` - High-dimensional neural network potential
* :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>` - registry-dependent interlayer potential (ILP)
* :doc:`ilp/tmd <pair_ilp_tmd>` - interlayer potential (ILP) potential for transition metal dichalcogenides (TMD)
* :doc:`kim <pair_kim>` - interface to potentials provided by KIM project
* :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>` - Kolmogorov-Crespi (KC) potential with no simplifications
* :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>` - Kolmogorov-Crespi (KC) potential with normals along z-axis
@ -305,6 +306,7 @@ accelerated styles exist.
* :doc:`reaxff <pair_reaxff>` - ReaxFF potential
* :doc:`rebo <pair_airebo>` - second generation REBO potential of Brenner
* :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
* :doc:`smd/hertz <pair_smd_hertz>` -
* :doc:`smd/tlsph <pair_smd_tlsph>` -

View File

@ -29,8 +29,8 @@ Syntax
.. parsed-literal::
*maxdelcs* value = delta1 delta2 (optional)
delta1 = The minimum thershold for cosine of three-body angle
delta2 = The maximum threshold for cosine of three-body angle
delta1 = The minimum thershold for the variation of cosine of three-body angle
delta2 = The maximum threshold for the variation of cosine of three-body angle
Examples
""""""""
@ -69,7 +69,7 @@ and K of atom I within a cutoff distance :math:`a `\sigma`.
The *sw/mod* style is designed for simulations of materials when
distinguishing three-body angles are necessary, such as borophene
and transition metal dichalcogenide, which cannot be described
and transition metal dichalcogenides, which cannot be described
by the original code for the Stillinger-Weber potential.
For instance, there are several types of angles around each Mo atom in `MoS_2`,
and some unnecessary angle types should be excluded in the three-body interaction.
@ -99,7 +99,7 @@ This smoothly turns off the energy and force contributions for :math:`\left| \de
It is suggested that :math:`\delta 1` and :math:`\delta_2` to be the value around
:math:`0.5 \left| \cos \theta_1 - \cos \theta_2 \right|`, with
:math:`\theta_1` and :math:`\theta_2` as the different types of angles around an atom.
For borophene and transition metal dichalcogenide, :math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`.
For borophene and transition metal dichalcogenides, :math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`.
This value enables the cut-off function to exclude unnecessary angles in the three-body SW terms.
.. note::

View File

@ -692,7 +692,7 @@ diagonalizers
diagonalizing
Diallo
diblock
dichalcogenide
dichalcogenides
Dickel
diel
Dietz
@ -2954,6 +2954,7 @@ safezone
Safran
Sagui
Saidi
saip
Salanne
Salles
sandia

View File

@ -0,0 +1 @@
../../../../potentials/BNCH.ILP

View File

@ -0,0 +1 @@
../../../../potentials/CH.rebo

View File

@ -0,0 +1 @@
../../../../potentials/MoS2.ILP

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,36 @@
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./bilayer_MoS2_AAprime_stacking.data
mass * 32.065 # mass of sulphur atom , uint: a.u.=1.66X10^(-27)kg
mass 1 95.94 # mass of molebdenum atom , uint: a.u.=1.66X10^(-27)kg
mass 4 95.94
# Define potentials
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
# Calculate the pair potential
compute 0 all pair ilp/tmd
compute 1 all pair sw/mod 1
compute 2 all pair sw/mod 2
variable SW1 equal c_1
variable SW2 equal c_2
variable ILP equal c_0
variable Eatt equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step etotal pe ke v_SW1 v_SW2 v_ILP v_Erep v_Eatt temp
thermo 100
thermo_modify lost error
timestep 1e-3
velocity all create 300.0 12345
fix intall all nve
run 1000

View File

@ -0,0 +1,154 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./bilayer_MoS2_AAprime_stacking.data
Reading data file ...
triclinic box = (0 0 -100) to (51.15232 44.299209 100) with tilt (25.57616 0 0)
1 by 1 by 1 MPI processor grid
reading atoms ...
1536 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.007 seconds
mass * 32.065 # mass of sulphur atom , uint: a.u.=1.66X10^(-27)kg
mass 1 95.94 # mass of molebdenum atom , uint: a.u.=1.66X10^(-27)kg
mass 4 95.94
# Define potentials
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
Reading ilp/tmd potential file MoS2.ILP with DATE: 2021-12-02
# Calculate the pair potential
compute 0 all pair ilp/tmd
compute 1 all pair sw/mod 1
compute 2 all pair sw/mod 2
variable SW1 equal c_1
variable SW2 equal c_2
variable ILP equal c_0
variable Eatt equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step etotal pe ke v_SW1 v_SW2 v_ILP v_Erep v_Eatt temp
thermo 100
thermo_modify lost error
timestep 1e-3
velocity all create 300.0 12345
fix intall all nve
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@Article{Ouyang2018
author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},
title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},
journal = {Nano Letters},
volume = 18,
pages = {6009}
year = 2018,
}
- ilp/tmd potential doi/10.1021/acs.jctc.1c00782
@Article{Ouyang2021
author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},
title = {Anisotropic Interlayer Force Field for Transition Metal Dichalcogenides: The Case of Molybdenum Disulfide},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {72377245}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 9 5 23
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(2) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(3) pair ilp/tmd, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual, copy from (3)
attributes: full, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 30.29 | 30.29 | 30.29 Mbytes
Step TotEng PotEng KinEng v_SW1 v_SW2 v_ILP v_Erep v_Eatt Temp
0 -1834.4469 -1893.9712 59.524297 -929.02881 -929.02881 -35.913549 63.00343 -98.916979 300
100 -1834.4497 -1883.3105 48.860775 -924.84264 -925.08505 -33.382796 56.58255 -89.965346 246.25629
200 -1834.4381 -1873.7072 39.269085 -922.29961 -922.52535 -28.882252 50.08277 -78.965022 197.91457
300 -1834.4483 -1881.1263 46.678028 -923.39264 -923.65627 -34.077402 51.011967 -85.089369 235.25534
400 -1834.431 -1868.0728 33.64182 -916.85743 -916.27044 -34.944916 50.414038 -85.358954 169.55338
500 -1834.4499 -1881.9059 47.456 -925.22919 -924.29582 -32.380877 44.755168 -77.136045 239.17628
600 -1834.4343 -1869.8642 35.429976 -920.97805 -919.60784 -29.278358 43.270241 -72.548599 178.56562
700 -1834.443 -1878.144 43.700934 -921.8051 -921.55671 -34.782141 49.959943 -84.742084 220.2509
800 -1834.4327 -1869.824 35.391298 -917.19193 -917.91383 -34.718247 55.349728 -90.067976 178.37068
900 -1834.4465 -1877.6741 43.227638 -923.33877 -922.50874 -31.826599 53.965592 -85.792191 217.86551
1000 -1834.4412 -1876.1808 41.739587 -923.17282 -923.49367 -29.514347 55.454643 -84.96899 210.3658
Loop time of 72.8218 on 1 procs for 1000 steps with 1536 atoms
Performance: 1.186 ns/day, 20.228 hours/ns, 13.732 timesteps/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 | 72.781 | 72.781 | 72.781 | 0.0 | 99.94
Bond | 0.00014503 | 0.00014503 | 0.00014503 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.016691 | 0.016691 | 0.016691 | 0.0 | 0.02
Output | 0.00057989 | 0.00057989 | 0.00057989 | 0.0 | 0.00
Modify | 0.013405 | 0.013405 | 0.013405 | 0.0 | 0.02
Other | | 0.01044 | | | 0.01
Nlocal: 1536 ave 1536 max 1536 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3510 ave 3510 max 3510 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 992256 ave 992256 max 992256 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 992256
Ave neighs/atom = 646
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:01:12

View File

@ -0,0 +1,154 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./bilayer_MoS2_AAprime_stacking.data
Reading data file ...
triclinic box = (0 0 -100) to (51.15232 44.299209 100) with tilt (25.57616 0 0)
2 by 2 by 1 MPI processor grid
reading atoms ...
1536 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.007 seconds
mass * 32.065 # mass of sulphur atom , uint: a.u.=1.66X10^(-27)kg
mass 1 95.94 # mass of molebdenum atom , uint: a.u.=1.66X10^(-27)kg
mass 4 95.94
# Define potentials
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
Reading sw potential file tmd.sw.mod with DATE: 2018-03-26
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
Reading ilp/tmd potential file MoS2.ILP with DATE: 2021-12-02
# Calculate the pair potential
compute 0 all pair ilp/tmd
compute 1 all pair sw/mod 1
compute 2 all pair sw/mod 2
variable SW1 equal c_1
variable SW2 equal c_2
variable ILP equal c_0
variable Eatt equal c_0[1]
variable Erep equal c_0[2]
thermo_style custom step etotal pe ke v_SW1 v_SW2 v_ILP v_Erep v_Eatt temp
thermo 100
thermo_modify lost error
timestep 1e-3
velocity all create 300.0 12345
fix intall all nve
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@Article{Ouyang2018
author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},
title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},
journal = {Nano Letters},
volume = 18,
pages = {6009}
year = 2018,
}
- ilp/tmd potential doi/10.1021/acs.jctc.1c00782
@Article{Ouyang2021
author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},
title = {Anisotropic Interlayer Force Field for Transition Metal Dichalcogenides: The Case of Molybdenum Disulfide},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {72377245}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 9 5 23
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(2) pair sw/mod, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(3) pair ilp/tmd, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual, copy from (3)
attributes: full, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 17.98 | 17.98 | 17.98 Mbytes
Step TotEng PotEng KinEng v_SW1 v_SW2 v_ILP v_Erep v_Eatt Temp
0 -1834.4469 -1893.9712 59.524297 -929.02881 -929.02881 -35.913549 63.00343 -98.916979 300
100 -1834.4497 -1883.3105 48.860775 -924.84264 -925.08505 -33.382796 56.58255 -89.965346 246.25629
200 -1834.4381 -1873.7072 39.269085 -922.29961 -922.52535 -28.882252 50.08277 -78.965022 197.91457
300 -1834.4483 -1881.1263 46.678028 -923.39264 -923.65627 -34.077402 51.011967 -85.089369 235.25534
400 -1834.431 -1868.0728 33.64182 -916.85743 -916.27044 -34.944916 50.414038 -85.358954 169.55338
500 -1834.4499 -1881.9059 47.456 -925.22919 -924.29582 -32.380877 44.755168 -77.136045 239.17628
600 -1834.4343 -1869.8642 35.429976 -920.97805 -919.60784 -29.278358 43.270241 -72.548599 178.56562
700 -1834.443 -1878.144 43.700934 -921.8051 -921.55671 -34.782141 49.959943 -84.742084 220.2509
800 -1834.4327 -1869.824 35.391298 -917.19193 -917.91383 -34.718247 55.349728 -90.067976 178.37068
900 -1834.4465 -1877.6741 43.227638 -923.33877 -922.50874 -31.826599 53.965592 -85.792191 217.86551
1000 -1834.4412 -1876.1808 41.739587 -923.17282 -923.49367 -29.514347 55.454643 -84.96899 210.3658
Loop time of 24.6309 on 4 procs for 1000 steps with 1536 atoms
Performance: 3.508 ns/day, 6.842 hours/ns, 40.599 timesteps/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 22.906 | 23.627 | 24.36 | 14.0 | 95.92
Bond | 0.00010889 | 0.00011572 | 0.00012719 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.25886 | 0.992 | 1.7126 | 68.3 | 4.03
Output | 0.00050901 | 0.00054202 | 0.00062029 | 0.0 | 0.00
Modify | 0.0030735 | 0.0033236 | 0.0035581 | 0.3 | 0.01
Other | | 0.008194 | | | 0.03
Nlocal: 384 ave 387 max 381 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Nghost: 2262 ave 2265 max 2259 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 248064 ave 250002 max 246126 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Total # of neighbors = 992256
Ave neighs/atom = 646
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:24

View File

@ -0,0 +1,143 @@
# DATE: 2018-03-26 UNITS: metal CONTRIBUTOR: Jin-Wu Jiang jwjiang5918@hotmail.com
# CITATION: J.-W. Jiang, Acta Mech. Solida. Sin 32, 17 (2019).
# The Stillinger-Weber parameters, for transition-metal dichalcogenide (TMD) lateral heterostructures.
# M = Mo, W; X = S, Se, Te
# these entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
# format of a single entry (one or more lines):
# element 1, element 2, element 3,
# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
# M-X-X terms
Mo S S 1.000 1.252 2.523 67.883 1.000 0.143 6.918 7.223 4 0 0.0
Mo Se Se 1.000 0.913 3.672 32.526 1.000 0.143 5.737 27.084 4 0 0.0
Mo Te Te 1.000 0.880 4.097 23.705 1.000 0.143 5.086 40.810 4 0 0.0
W S S 1.000 0.889 3.558 37.687 1.000 0.143 5.664 24.525 4 0 0.0
W Se Se 1.000 0.706 4.689 25.607 1.000 0.143 5.476 65.662 4 0 0.0
W Te Te 1.000 0.778 4.632 21.313 1.000 0.143 4.326 62.148 4 0 0.0
# X-M-M terms
S Mo Mo 1.000 1.252 2.523 62.449 1.000 0.143 6.918 7.223 4 0 0.0
S W W 1.000 0.889 3.558 33.553 1.000 0.143 5.664 24.525 4 0 0.0
Se Mo Mo 1.000 0.913 3.672 27.079 1.000 0.143 5.737 27.084 4 0 0.0
Se W W 1.000 0.706 4.689 23.218 1.000 0.143 5.476 65.662 4 0 0.0
Te Mo Mo 1.000 0.880 4.097 20.029 1.000 0.143 5.086 40.810 4 0 0.0
Te W W 1.000 0.778 4.632 17.370 1.000 0.143 4.326 62.148 4 0 0.0
# M-X1-X2 terms
Mo S Se 1.000 0.000 0.000 46.989 1.000 0.143 0.000 0.000 4 0 0.0
Mo S Te 1.000 0.000 0.000 40.114 1.000 0.143 0.000 0.000 4 0 0.0
Mo Se S 1.000 0.000 0.000 46.989 1.000 0.143 0.000 0.000 4 0 0.0
Mo Se Te 1.000 0.000 0.000 27.767 1.000 0.143 0.000 0.000 4 0 0.0
Mo Te S 1.000 0.000 0.000 40.114 1.000 0.143 0.000 0.000 4 0 0.0
Mo Te Se 1.000 0.000 0.000 27.767 1.000 0.143 0.000 0.000 4 0 0.0
W S Se 1.000 0.000 0.000 31.065 1.000 0.143 0.000 0.000 4 0 0.0
W S Te 1.000 0.000 0.000 28.341 1.000 0.143 0.000 0.000 4 0 0.0
W Se S 1.000 0.000 0.000 31.065 1.000 0.143 0.000 0.000 4 0 0.0
W Se Te 1.000 0.000 0.000 23.362 1.000 0.143 0.000 0.000 4 0 0.0
W Te S 1.000 0.000 0.000 28.341 1.000 0.143 0.000 0.000 4 0 0.0
W Te Se 1.000 0.000 0.000 23.362 1.000 0.143 0.000 0.000 4 0 0.0
# X-M1-M2 terms
S Mo W 1.000 0.000 0.000 45.775 1.000 0.143 0.000 0.000 4 0 0.0
S W Mo 1.000 0.000 0.000 45.775 1.000 0.143 0.000 0.000 4 0 0.0
Se Mo W 1.000 0.000 0.000 25.074 1.000 0.143 0.000 0.000 4 0 0.0
Se W Mo 1.000 0.000 0.000 25.074 1.000 0.143 0.000 0.000 4 0 0.0
Te Mo W 1.000 0.000 0.000 18.652 1.000 0.143 0.000 0.000 4 0 0.0
Te W Mo 1.000 0.000 0.000 18.652 1.000 0.143 0.000 0.000 4 0 0.0
# zero terms
Mo Mo Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Mo Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
W Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S S Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Se Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
S Te Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se S Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Se Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Se Te Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Mo S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Mo Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Mo Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te W S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te W Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te W Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te S Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Se Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te Mo 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te W 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te S 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te Se 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0
Te Te Te 0.000 1.000 1.000 1.000 1.000 0.143 1.000 1.000 4 0 0.0

View File

@ -0,0 +1,218 @@
Gold/graphene heterjunction atop2
206 atoms
2 atom types
0 17.216640 xlo xhi
0 14.910048 ylo yhi
-30 30.000000 zlo zhi
8.60832 0 0 xy xz yz
Atoms
1 1 1 0 -8.60832 -4.97002 3.5
2 1 1 0 -7.1736 -4.14166 5.85461
3 1 1 0 -5.73888 -3.31334 8.20922
4 1 1 0 -7.1736 -2.48501 3.5
5 1 1 0 -5.73888 -1.65666 5.85461
6 1 1 0 -4.30416 -0.828336 8.20922
7 1 1 0 -5.73888 0 3.5
8 1 1 0 -4.30416 0.828352 5.85461
9 1 1 0 -2.86944 1.65667 8.20922
10 1 1 0 -4.30416 2.48501 3.5
11 1 1 0 -2.86944 3.31336 5.85461
12 1 1 0 -1.43472 4.14168 8.20922
13 1 1 0 -2.86944 4.97002 3.5
14 1 1 0 -1.43472 5.79837 5.85461
15 1 1 0 0 6.62669 8.20922
16 1 1 0 -1.43472 7.45502 3.5
17 1 1 0 0 8.28337 5.85461
18 1 1 0 1.43472 9.11167 8.20922
19 1 1 0 -5.73888 -4.97002 3.5
20 1 1 0 -4.30416 -4.14166 5.85461
21 1 1 0 -2.86944 -3.31334 8.20922
22 1 1 0 -4.30416 -2.48501 3.5
23 1 1 0 -2.86944 -1.65666 5.85461
24 1 1 0 -1.43472 -0.828336 8.20922
25 1 1 0 -2.86944 0 3.5
26 1 1 0 -1.43472 0.828352 5.85461
27 1 1 0 0 1.65667 8.20922
28 1 1 0 -1.43472 2.48501 3.5
29 1 1 0 0 3.31336 5.85461
30 1 1 0 1.43472 4.14168 8.20922
31 1 1 0 0 4.97002 3.5
32 1 1 0 1.43472 5.79837 5.85461
33 1 1 0 2.86944 6.62669 8.20922
34 1 1 0 1.43472 7.45502 3.5
35 1 1 0 2.86944 8.28337 5.85461
36 1 1 0 4.30416 9.11167 8.20922
37 1 1 0 -2.86944 -4.97002 3.5
38 1 1 0 -1.43472 -4.14166 5.85461
39 1 1 0 0 -3.31334 8.20922
40 1 1 0 -1.43472 -2.48501 3.5
41 1 1 0 0 -1.65666 5.85461
42 1 1 0 1.43472 -0.828336 8.20922
43 1 1 0 0 0 3.5
44 1 1 0 1.43472 0.828352 5.85461
45 1 1 0 2.86944 1.65667 8.20922
46 1 1 0 1.43472 2.48501 3.5
47 1 1 0 2.86944 3.31336 5.85461
48 1 1 0 4.30416 4.14168 8.20922
49 1 1 0 2.86944 4.97002 3.5
50 1 1 0 4.30416 5.79837 5.85461
51 1 1 0 5.73888 6.62669 8.20922
52 1 1 0 4.30416 7.45502 3.5
53 1 1 0 5.73888 8.28337 5.85461
54 1 1 0 7.1736 9.11167 8.20922
55 1 1 0 0 -4.97002 3.5
56 1 1 0 1.43472 -4.14166 5.85461
57 1 1 0 2.86944 -3.31334 8.20922
58 1 1 0 1.43472 -2.48501 3.5
59 1 1 0 2.86944 -1.65666 5.85461
60 1 1 0 4.30416 -0.828336 8.20922
61 1 1 0 2.86944 0 3.5
62 1 1 0 4.30416 0.828352 5.85461
63 1 1 0 5.73888 1.65667 8.20922
64 1 1 0 4.30416 2.48501 3.5
65 1 1 0 5.73888 3.31336 5.85461
66 1 1 0 7.1736 4.14168 8.20922
67 1 1 0 5.73888 4.97002 3.5
68 1 1 0 7.1736 5.79837 5.85461
69 1 1 0 8.60832 6.62669 8.20922
70 1 1 0 7.1736 7.45502 3.5
71 1 1 0 8.60832 8.28337 5.85461
72 1 1 0 10.0431 9.11167 8.20922
73 1 1 0 2.86944 -4.97002 3.5
74 1 1 0 4.30416 -4.14166 5.85461
75 1 1 0 5.73888 -3.31334 8.20922
76 1 1 0 4.30416 -2.48501 3.5
77 1 1 0 5.73888 -1.65666 5.85461
78 1 1 0 7.1736 -0.828336 8.20922
79 1 1 0 5.73888 0 3.5
80 1 1 0 7.1736 0.828352 5.85461
81 1 1 0 8.60832 1.65667 8.20922
82 1 1 0 7.1736 2.48501 3.5
83 1 1 0 8.60832 3.31336 5.85461
84 1 1 0 10.0431 4.14168 8.20922
85 1 1 0 8.60832 4.97002 3.5
86 1 1 0 10.0431 5.79837 5.85461
87 1 1 0 11.4777 6.62669 8.20922
88 1 1 0 10.0431 7.45502 3.5
89 1 1 0 11.4777 8.28337 5.85461
90 1 1 0 12.9124 9.11167 8.20922
91 1 1 0 5.73888 -4.97002 3.5
92 1 1 0 7.1736 -4.14166 5.85461
93 1 1 0 8.60832 -3.31334 8.20922
94 1 1 0 7.1736 -2.48501 3.5
95 1 1 0 8.60832 -1.65666 5.85461
96 1 1 0 10.0431 -0.828336 8.20922
97 1 1 0 8.60832 0 3.5
98 1 1 0 10.0431 0.828352 5.85461
99 1 1 0 11.4777 1.65667 8.20922
100 1 1 0 10.0431 2.48501 3.5
101 1 1 0 11.4777 3.31336 5.85461
102 1 1 0 12.9124 4.14168 8.20922
103 1 1 0 11.4777 4.97002 3.5
104 1 1 0 12.9124 5.79837 5.85461
105 1 1 0 14.3472 6.62669 8.20922
106 1 1 0 12.9124 7.45502 3.5
107 1 1 0 14.3472 8.28337 5.85461
108 1 1 0 15.7819 9.11167 8.20922
109 2 2 0 -11.0678 -6.39002 0
110 2 2 0 -9.83808 -5.68001 0
111 2 2 0 -9.83808 -4.26001 0
112 2 2 0 -8.60832 -3.55001 0
113 2 2 0 -8.60832 -2.13001 0
114 2 2 0 -7.37856 -1.42 0
115 2 2 0 -7.37856 0 0
116 2 2 0 -6.1488 0.710007 0
117 2 2 0 -6.1488 2.13001 0
118 2 2 0 -4.91904 2.84001 0
119 2 2 0 -4.91904 4.26001 0
120 2 2 0 -3.68928 4.97002 0
121 2 2 0 -3.68928 6.39002 0
122 2 2 0 -2.45952 7.10003 0
123 2 2 0 -8.60832 -6.39002 0
124 2 2 0 -7.37856 -5.68001 0
125 2 2 0 -7.37856 -4.26001 0
126 2 2 0 -6.1488 -3.55001 0
127 2 2 0 -6.1488 -2.13001 0
128 2 2 0 -4.91904 -1.42 0
129 2 2 0 -4.91904 0 0
130 2 2 0 -3.68928 0.710007 0
131 2 2 0 -3.68928 2.13001 0
132 2 2 0 -2.45952 2.84001 0
133 2 2 0 -2.45952 4.26001 0
134 2 2 0 -1.22976 4.97002 0
135 2 2 0 -1.22976 6.39002 0
136 2 2 0 0 7.10003 0
137 2 2 0 -6.1488 -6.39002 0
138 2 2 0 -4.91904 -5.68001 0
139 2 2 0 -4.91904 -4.26001 0
140 2 2 0 -3.68928 -3.55001 0
141 2 2 0 -3.68928 -2.13001 0
142 2 2 0 -2.45952 -1.42 0
143 2 2 0 -2.45952 0 0
144 2 2 0 -1.22976 0.710007 0
145 2 2 0 -1.22976 2.13001 0
146 2 2 0 0 2.84001 0
147 2 2 0 0 4.26001 0
148 2 2 0 1.22976 4.97002 0
149 2 2 0 1.22976 6.39002 0
150 2 2 0 2.45952 7.10003 0
151 2 2 0 -3.68928 -6.39002 0
152 2 2 0 -2.45952 -5.68001 0
153 2 2 0 -2.45952 -4.26001 0
154 2 2 0 -1.22976 -3.55001 0
155 2 2 0 -1.22976 -2.13001 0
156 2 2 0 0 -1.42 0
157 2 2 0 0 0 0
158 2 2 0 1.22976 0.710007 0
159 2 2 0 1.22976 2.13001 0
160 2 2 0 2.45952 2.84001 0
161 2 2 0 2.45952 4.26001 0
162 2 2 0 3.68928 4.97002 0
163 2 2 0 3.68928 6.39002 0
164 2 2 0 4.91904 7.10003 0
165 2 2 0 -1.22976 -6.39002 0
166 2 2 0 0 -5.68001 0
167 2 2 0 0 -4.26001 0
168 2 2 0 1.22976 -3.55001 0
169 2 2 0 1.22976 -2.13001 0
170 2 2 0 2.45952 -1.42 0
171 2 2 0 2.45952 0 0
172 2 2 0 3.68928 0.710007 0
173 2 2 0 3.68928 2.13001 0
174 2 2 0 4.91904 2.84001 0
175 2 2 0 4.91904 4.26001 0
176 2 2 0 6.1488 4.97002 0
177 2 2 0 6.1488 6.39002 0
178 2 2 0 7.37856 7.10003 0
179 2 2 0 1.22976 -6.39002 0
180 2 2 0 2.45952 -5.68001 0
181 2 2 0 2.45952 -4.26001 0
182 2 2 0 3.68928 -3.55001 0
183 2 2 0 3.68928 -2.13001 0
184 2 2 0 4.91904 -1.42 0
185 2 2 0 4.91904 0 0
186 2 2 0 6.1488 0.710007 0
187 2 2 0 6.1488 2.13001 0
188 2 2 0 7.37856 2.84001 0
189 2 2 0 7.37856 4.26001 0
190 2 2 0 8.60832 4.97002 0
191 2 2 0 8.60832 6.39002 0
192 2 2 0 9.83808 7.10003 0
193 2 2 0 3.68928 -6.39002 0
194 2 2 0 4.91904 -5.68001 0
195 2 2 0 4.91904 -4.26001 0
196 2 2 0 6.1488 -3.55001 0
197 2 2 0 6.1488 -2.13001 0
198 2 2 0 7.37856 -1.42 0
199 2 2 0 7.37856 0 0
200 2 2 0 8.60832 0.710007 0
201 2 2 0 8.60832 2.13001 0
202 2 2 0 9.83808 2.84001 0
203 2 2 0 9.83808 4.26001 0
204 2 2 0 11.0678 4.97002 0
205 2 2 0 11.0678 6.39002 0
206 2 2 0 12.2976 7.10003 0

View File

@ -0,0 +1 @@
../../../../potentials/Au_u3.eam

View File

@ -0,0 +1 @@
../../../../potentials/CH.rebo

View File

@ -0,0 +1 @@
../../../../potentials/CHAu.ILP

View File

@ -0,0 +1,44 @@
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./3Lgold_1Lgr_atop_sliding.data
# global group definition
group gold type 1
group gra type 2
# Define mass
mass * 12.0107 # mass of carbon atom , uint: a.u.=1.66X10^(-27)kg
mass 1 196.96657 # mass of gold atom , uint: a.u.=1.66X10^(-27)kg
# Define potentials
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam ./Au_u3.eam
pair_coeff * * rebo ./CH.rebo NULL C
pair_coeff * * saip/metal ./CHAu.ILP Au C
# compute energy
compute 0 all pair rebo
compute 1 all pair eam
compute 2 all pair saip/metal
variable REBO equal c_0
variable EAM equal c_1
variable ILP equal c_2
thermo_style custom step etotal pe ke v_REBO v_ILP temp
thermo 100
thermo_modify lost error
# Creat initial velocity
velocity all set 0.0 0.0 0.0
velocity gra create 300.0 4928459 mom yes rot yes dist gaussian
velocity gold create 300.0 4928459 mom yes rot yes dist gaussian
# Integration
fix intsub gold nve
fix intrib gra nve
timestep 1e-3
run 1000

View File

@ -0,0 +1,164 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./3Lgold_1Lgr_atop_sliding.data
Reading data file ...
triclinic box = (0 0 -30) to (17.21664 14.910048 30) with tilt (8.60832 0 0)
1 by 1 by 1 MPI processor grid
reading atoms ...
206 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.003 seconds
# global group definition
group gold type 1
108 atoms in group gold
group gra type 2
98 atoms in group gra
# Define mass
mass * 12.0107 # mass of carbon atom , uint: a.u.=1.66X10^(-27)kg
mass 1 196.96657 # mass of gold atom , uint: a.u.=1.66X10^(-27)kg
# Define potentials
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam ./Au_u3.eam
Reading eam potential file ./Au_u3.eam with DATE: 2007-06-11
pair_coeff * * rebo ./CH.rebo NULL C
Reading rebo potential file ./CH.rebo with DATE: 2018-7-3
pair_coeff * * saip/metal ./CHAu.ILP Au C
Reading saip/metal potential file ./CHAu.ILP with DATE: 2021-12-02
# compute energy
compute 0 all pair rebo
compute 1 all pair eam
compute 2 all pair saip/metal
variable REBO equal c_0
variable EAM equal c_1
variable ILP equal c_2
thermo_style custom step etotal pe ke v_REBO v_ILP temp
thermo 100
thermo_modify lost error
# Creat initial velocity
velocity all set 0.0 0.0 0.0
velocity gra create 300.0 4928459 mom yes rot yes dist gaussian
velocity gold create 300.0 4928459 mom yes rot yes dist gaussian
# Integration
fix intsub gold nve
fix intrib gra nve
timestep 1e-3
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@Article{Ouyang2018
author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},
title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},
journal = {Nano Letters},
volume = 18,
pages = {6009}
year = 2018,
}
- saip/metal potential doi.org/10.1021/acs.jctc.1c00622
@Article{Ouyang2021
author = {W. Ouyang, O. Hod, and R. Guerra},
title = {Registry-Dependent Potential for Interfaces of Gold with Graphitic Systems},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {7215-7223}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 3 2 7
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair eam, perpetual, skip from (4)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(2) pair rebo, perpetual, skip from (3)
attributes: full, newton on, ghost
pair build: skip/ghost
stencil: none
bin: none
(3) pair saip/metal, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 9.747 | 9.747 | 9.747 Mbytes
Step TotEng PotEng KinEng v_REBO v_ILP Temp
0 -1121.4621 -1129.3728 7.9107209 -724.70925 -6.0302289 298.53659
100 -1121.4483 -1124.9731 3.5248501 -723.03272 -5.9765533 133.02159
200 -1121.4403 -1125.2912 3.8509646 -722.66784 -6.0468507 145.32858
300 -1121.4424 -1126.4665 5.0240531 -722.72787 -6.0350568 189.59886
400 -1121.4419 -1125.1443 3.7023978 -722.59976 -5.8294548 139.72193
500 -1121.4413 -1125.2711 3.8297939 -722.5342 -6.0189944 144.52963
600 -1121.4449 -1125.8808 4.4359049 -722.77707 -5.8685221 167.40319
700 -1121.4489 -1126.1966 4.747709 -723.13681 -5.8666379 179.17012
800 -1121.4443 -1125.8469 4.402607 -722.94487 -6.0094567 166.14658
900 -1121.4424 -1125.8437 4.4013317 -722.94918 -5.8699702 166.09846
1000 -1121.4467 -1125.7453 4.2986881 -722.66682 -6.0651168 162.22487
Loop time of 6.43246 on 1 procs for 1000 steps with 206 atoms
Performance: 13.432 ns/day, 1.787 hours/ns, 155.462 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 6.4201 | 6.4201 | 6.4201 | 0.0 | 99.81
Bond | 8.9059e-05 | 8.9059e-05 | 8.9059e-05 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0071852 | 0.0071852 | 0.0071852 | 0.0 | 0.11
Output | 0.00026031 | 0.00026031 | 0.00026031 | 0.0 | 0.00
Modify | 0.0019433 | 0.0019433 | 0.0019433 | 0.0 | 0.03
Other | | 0.002875 | | | 0.04
Nlocal: 206 ave 206 max 206 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 2187 ave 2187 max 2187 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 158548 ave 158548 max 158548 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 158548
Ave neighs/atom = 769.65049
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:06

View File

@ -0,0 +1,164 @@
LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
units metal
atom_style full
processors * * 1
boundary p p f
read_data ./3Lgold_1Lgr_atop_sliding.data
Reading data file ...
triclinic box = (0 0 -30) to (17.21664 14.910048 30) with tilt (8.60832 0 0)
2 by 2 by 1 MPI processor grid
reading atoms ...
206 atoms
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.003 seconds
# global group definition
group gold type 1
108 atoms in group gold
group gra type 2
98 atoms in group gra
# Define mass
mass * 12.0107 # mass of carbon atom , uint: a.u.=1.66X10^(-27)kg
mass 1 196.96657 # mass of gold atom , uint: a.u.=1.66X10^(-27)kg
# Define potentials
pair_style hybrid/overlay eam rebo saip/metal 16.0
pair_coeff 1 1 eam ./Au_u3.eam
Reading eam potential file ./Au_u3.eam with DATE: 2007-06-11
pair_coeff * * rebo ./CH.rebo NULL C
Reading rebo potential file ./CH.rebo with DATE: 2018-7-3
pair_coeff * * saip/metal ./CHAu.ILP Au C
Reading saip/metal potential file ./CHAu.ILP with DATE: 2021-12-02
# compute energy
compute 0 all pair rebo
compute 1 all pair eam
compute 2 all pair saip/metal
variable REBO equal c_0
variable EAM equal c_1
variable ILP equal c_2
thermo_style custom step etotal pe ke v_REBO v_ILP temp
thermo 100
thermo_modify lost error
# Creat initial velocity
velocity all set 0.0 0.0 0.0
velocity gra create 300.0 4928459 mom yes rot yes dist gaussian
velocity gold create 300.0 4928459 mom yes rot yes dist gaussian
# Integration
fix intsub gold nve
fix intrib gra nve
timestep 1e-3
run 1000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848
@Article{Ouyang2018
author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},
title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},
journal = {Nano Letters},
volume = 18,
pages = {6009}
year = 2018,
}
- saip/metal potential doi.org/10.1021/acs.jctc.1c00622
@Article{Ouyang2021
author = {W. Ouyang, O. Hod, and R. Guerra},
title = {Registry-Dependent Potential for Interfaces of Gold with Graphitic Systems},
journal = {J. Chem. Theory Comput.},
volume = 17,
pages = {7215-7223}
year = 2021,
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 18
ghost atom cutoff = 18
binsize = 9, bins = 3 2 7
4 neighbor lists, perpetual/occasional/extra = 4 0 0
(1) pair eam, perpetual, skip from (4)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
(2) pair rebo, perpetual, skip from (3)
attributes: full, newton on, ghost
pair build: skip/ghost
stencil: none
bin: none
(3) pair saip/metal, perpetual
attributes: full, newton on, ghost
pair build: full/bin/ghost
stencil: full/ghost/bin/3d
bin: standard
(4) neighbor class addition, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 9.299 | 9.299 | 9.3 Mbytes
Step TotEng PotEng KinEng v_REBO v_ILP Temp
0 -1121.4621 -1129.3728 7.9107209 -724.70925 -6.0302289 298.53659
100 -1121.4483 -1124.9731 3.5248501 -723.03272 -5.9765533 133.02159
200 -1121.4403 -1125.2912 3.8509646 -722.66784 -6.0468507 145.32858
300 -1121.4424 -1126.4665 5.0240531 -722.72787 -6.0350568 189.59886
400 -1121.4419 -1125.1443 3.7023978 -722.59976 -5.8294548 139.72193
500 -1121.4413 -1125.2711 3.8297939 -722.5342 -6.0189944 144.52963
600 -1121.4449 -1125.8808 4.4359049 -722.77707 -5.8685221 167.40319
700 -1121.4489 -1126.1966 4.747709 -723.13681 -5.8666379 179.17012
800 -1121.4443 -1125.8469 4.402607 -722.94487 -6.0094567 166.14658
900 -1121.4424 -1125.8437 4.4013317 -722.94918 -5.8699702 166.09846
1000 -1121.4467 -1125.7453 4.2986881 -722.66682 -6.0651168 162.22487
Loop time of 2.44776 on 4 procs for 1000 steps with 206 atoms
Performance: 35.298 ns/day, 0.680 hours/ns, 408.537 timesteps/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.6906 | 2.0172 | 2.3939 | 19.2 | 82.41
Bond | 5.4278e-05 | 6.3818e-05 | 7.3153e-05 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.049363 | 0.42514 | 0.75161 | 41.7 | 17.37
Output | 0.00018596 | 0.00020656 | 0.00024754 | 0.0 | 0.01
Modify | 0.00063656 | 0.00074701 | 0.00089514 | 0.0 | 0.03
Other | | 0.004362 | | | 0.18
Nlocal: 51.5 ave 61 max 44 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Nghost: 1670.75 ave 1699 max 1641 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 39637 ave 47080 max 33737 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Total # of neighbors = 158548
Ave neighs/atom = 769.65049
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:02

18
potentials/CHAu.ILP Normal file
View File

@ -0,0 +1,18 @@
# DATE: 2021-12-02 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com
# CITATION: W. Ouyang, O. Hod, and R. Guerra, J. Chem. Theory Comput. 17, 7215 (2021).
# Semi-anisotropic interlayer Potential (SAIP) for graphene/gold, benzene/gold heterojunctions
# The parameters below are fitted against the PBE + D3 DFT reference data.
# The parameters for C-C, C-H and H-H are taken from Nano Letters 18, 6009-6016 (2018).
#
# --------------- Repulsion Potential ----------------++++++++++++++ Vdw Potential ++++++++++++++++*****
# beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut
C C 3.205843 7.511126 1.235334 1.528338E-5 37.530428 15.499947 0.7954443 3.681440 25.714535E3 1.0 2.0
H H 3.974540 6.53799 1.080633 0.6700556 0.8333833 15.022371 0.7490632 2.767223 1.6159581E3 1.0 1.2
C H 2.642950 12.91410 1.020257 0.9750012 25.340996 15.222927 0.8115998 3.887324 5.6874617E3 1.0 1.5
Au C 3.6913278482 13.5655648421 1.0175514400 0.0070964784 -0.0010368264 11.0586486772 1.0635582839 3.7552608806 81.5847131142 1000.0 1.0
Au H 3.7899616131 4.3065009639 10.6811825156 0.2250887843 -0.1116891520 18.6149213872 0.9833194192 3.3507558896 70.6865381780 1000.0 1.0
Au Au 3.6671967387 12.8109735143 1.0353581041 0.0000000000 0.0000000000 10.1628585345 1.0642897301 3.7372959779 0.0000000000 1000.0 1.0
# Symmetric part
C Au 3.6913278482 13.5655648421 1.0175514400 0.0070964784 -0.0010368264 11.0586486772 1.0635582839 3.7552608806 81.5847131142 1000.0 2.0
H Au 3.7899616131 4.3065009639 10.6811825156 0.2250887843 -0.1116891520 18.6149213872 0.9833194192 3.3507558896 70.6865381780 1000.0 1.5
H C 2.642950 12.91410 1.020257 0.9750012 25.340996 15.222927 0.8115998 3.887324 5.6874617E3 1.0 1.5

12
potentials/MoS2.ILP Normal file
View File

@ -0,0 +1,12 @@
# DATE: 2021-12-02 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com
# CITATION: W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).
# Interlayer Potential (ILP) for bilayer and bulk Molybdenum Disulfide.
# The parameters below are fitted against the HSE + MBD-NL DFT reference data.
#
# --------------- Repulsion Potential ----------------++++++++++++++ Vdw Potential ++++++++++++++++*****
# beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut
Mo Mo 5.5794503728 9.3776624643 2.0272224617 144.1517754265 97.9785701736 89.4375967588 2.0590305447 5.1220549022 491850.3161946840 1.0000 4.0
S S 3.1614021873 8.0932627454 1.9531402037 4.5867641760 118.0654664611 58.8094158385 0.2153665787 4.2996001250 148811.2434089213 1.0000 4.0
Mo S 3.6271521213 19.9713750996 7.5850312329 76.1019310047 3.3174955929 45.7203284638 0.9474703731 4.4104250318 150597.8577162951 1.0000 4.0
# symmetric part
S Mo 3.6271521213 19.9713750996 7.5850312329 76.1019310047 3.3174955929 45.7203284638 0.9474703731 4.4104250318 150597.8577162951 1.0000 4.0

View File

@ -55,7 +55,8 @@ NPairSkipIntel::~NPairSkipIntel() {
void NPairSkipIntel::copy_neighbor_info()
{
NPair::copy_neighbor_info();
if (_full_props) delete []_full_props;
// Only need to set _full_props once; npair object deleted for changes
if (_full_props) return;
_full_props = new int[neighbor->nrequest];
for (int i = 0; i < neighbor->nrequest; i++)
_full_props[i] = neighbor->requests[i]->full;

View File

@ -37,6 +37,7 @@
#include <cmath>
#include <cstring>
#include <map>
using namespace LAMMPS_NS;
using namespace InterLayer;
@ -56,9 +57,15 @@ static const char cite_ilp[] =
" year = 2018,\n"
"}\n\n";
// to indicate which potential style was used in outputs
static std::map<int, std::string> variant_map = {
{PairILPGrapheneHBN::ILP_GrhBN, "ilp/graphene/hbn"},
{PairILPGrapheneHBN::ILP_TMD, "ilp/tmd"},
{PairILPGrapheneHBN::SAIP_METAL, "saip/metal"}};
/* ---------------------------------------------------------------------- */
PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp)
PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp), variant(ILP_GrhBN)
{
restartinfo = 0;
one_coeff = 1;
@ -86,6 +93,14 @@ PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp)
dnormal = nullptr;
dnormdri = nullptr;
// for ilp/tmd
dnn = nullptr;
vect = nullptr;
pvet = nullptr;
dpvet1 = nullptr;
dpvet2 = nullptr;
dNave = nullptr;
// always compute energy offset
offset_flag = 1;
@ -104,6 +119,13 @@ PairILPGrapheneHBN::~PairILPGrapheneHBN()
memory->destroy(normal);
memory->destroy(dnormal);
memory->destroy(dnormdri);
// adds for ilp/tmd
memory->destroy(dnn);
memory->destroy(vect);
memory->destroy(pvet);
memory->destroy(dpvet1);
memory->destroy(dpvet2);
memory->destroy(dNave);
if (allocated) {
memory->destroy(setflag);
@ -194,7 +216,7 @@ void PairILPGrapheneHBN::read_file(char *filename)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, filename, "ilp/graphene/hbn", unit_convert_flag);
PotentialFileReader reader(lmp, filename, variant_map[variant], unit_convert_flag);
char *line;
// transparently convert units for supported conversions
@ -292,11 +314,15 @@ void PairILPGrapheneHBN::read_file(char *filename)
int n = -1;
for (int m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement) {
if (n >= 0) error->all(FLERR, "ILP potential file has duplicate entry");
if (n >= 0)
error->all(FLERR, "{} potential file {} has a duplicate entry", variant_map[variant],
filename);
n = m;
}
}
if (n < 0) error->all(FLERR, "Potential file is missing an entry");
if (n < 0)
error->all(FLERR, "{} potential file {} is missing an entry", variant_map[variant],
filename);
elem2param[i][j] = n;
cutILPsq[i][j] = params[n].rcut * params[n].rcut;
}
@ -310,9 +336,9 @@ void PairILPGrapheneHBN::read_file(char *filename)
void PairILPGrapheneHBN::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR, "Pair style ilp/graphene/hbn requires newton pair on");
error->all(FLERR, "Pair style {} requires newton pair on", variant_map[variant]);
if (!atom->molecule_flag)
error->all(FLERR, "Pair style ilp/graphene/hbn requires atom attribute molecule");
error->all(FLERR, "Pair style {} requires atom attribute molecule", variant_map[variant]);
// need a full neighbor list, including neighbors of ghosts

View File

@ -34,16 +34,16 @@ class PairILPGrapheneHBN : public Pair {
void coeff(int, char **) override;
double init_one(int, int) override;
void init_style() override;
void ILP_neigh();
void calc_normal();
void calc_FRep(int, int);
void calc_FvdW(int, int);
double single(int, int, int, int, double, double, double, double &) override;
static constexpr int NPARAMS_PER_LINE = 13;
enum { ILP_GrhBN, ILP_TMD, SAIP_METAL }; // for telling class variants apart in shared code
protected:
int me;
int variant;
int maxlocal; // size of numneigh, firstneigh arrays
int pgsize; // size of neighbor page
int oneatom; // max # of neighbors for one atom
@ -68,6 +68,18 @@ class PairILPGrapheneHBN : public Pair {
double ***dnormdri;
double ****dnormal;
// adds for ilp/tmd
int Nnei; // max # of nearest neighbors for one atom
double **dnn;
double **vect;
double **pvet;
double ***dpvet1;
double ***dpvet2;
double ***dNave;
virtual void ILP_neigh();
virtual void calc_normal();
virtual void calc_FRep(int, int);
void read_file(char *);
void allocate();
};

View File

@ -0,0 +1,885 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
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: Wengen Ouyang (Wuhan University)
e-mail: w.g.ouyang at gmail dot com
This is a full version of the potential described in
[Ouyang et al, J. Chem. Theory Comput. 17, 7237 (2021).]
------------------------------------------------------------------------- */
#include "pair_ilp_tmd.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "interlayer_taper.h"
#include "memory.h"
#include "my_page.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace InterLayer;
#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
static const char cite_ilp_tmd[] = "ilp/tmd potential doi/10.1021/acs.jctc.1c00782\n"
"@Article{Ouyang2021\n"
" author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. "
"Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},\n"
" title = {Anisotropic Interlayer Force Field for Transition "
"Metal Dichalcogenides: The Case of Molybdenum Disulfide},\n"
" journal = {J. Chem. Theory Comput.},\n"
" volume = 17,\n"
" pages = {72377245}\n"
" year = 2021,\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairILPTMD::PairILPTMD(LAMMPS *lmp) : PairILPGrapheneHBN(lmp)
{
variant = ILP_TMD;
single_enable = 0;
// for TMD, each atom have six neighbors
Nnei = 6;
if (lmp->citeme) lmp->citeme->add(cite_ilp_tmd);
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairILPTMD::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command");
if (!utils::strmatch(force->pair_style, "^hybrid/overlay"))
error->all(FLERR, "Pair style ilp/tmd must be used as sub-style with hybrid/overlay");
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp);
}
/* ----------------------------------------------------------------------
Repulsive forces and energy
------------------------------------------------------------------------- */
void PairILPTMD::calc_FRep(int eflag, int /* vflag */)
{
int i, j, ii, jj, inum, jnum, itype, jtype, k, kk;
double prodnorm1, fkcx, fkcy, fkcz;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1;
double rsq, r, Rcut, rhosq1, exp0, exp1, r2inv, r6inv, r8inv, Tap, dTap, Vilp;
double frho1, TSvdw, TSvdw2inv, Erep, fsum, rdsq1;
int *ilist, *jlist, *numneigh, **firstneigh;
int *ILP_neighs_i;
evdwl = 0.0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double dprodnorm1[3] = {0.0, 0.0, 0.0};
double fp1[3] = {0.0, 0.0, 0.0};
double fprod1[3] = {0.0, 0.0, 0.0};
double delki[3] = {0.0, 0.0, 0.0};
double fk[3] = {0.0, 0.0, 0.0};
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
//calculate exp(-lambda*(r-z0))*[epsilon/2 + f(rho_ij)]
// loop over neighbors of owned atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// only include the interation between different layers
if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
int iparam_ij = elem2param[map[itype]][map[jtype]];
Param &p = params[iparam_ij];
r = sqrt(rsq);
// turn on/off taper function
if (tap_flag) {
Rcut = sqrt(cutsq[itype][jtype]);
Tap = calc_Tap(r, Rcut);
dTap = calc_dTap(r, Rcut);
} else {
Tap = 1.0;
dTap = 0.0;
}
// Calculate the transverse distance
// note that rho_ij does not equal to rho_ji except when normals are all along z
prodnorm1 = normal[i][0] * delx + normal[i][1] * dely + normal[i][2] * delz;
rhosq1 = rsq - prodnorm1 * prodnorm1; // rho_ij
rdsq1 = rhosq1 * p.delta2inv; // (rho_ij/delta)^2
// store exponents
exp0 = exp(-p.lambda * (r - p.z0));
exp1 = exp(-rdsq1);
frho1 = exp1 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda * exp0 / r * Erep;
fpair1 = 2.0 * exp0 * frho1 * p.delta2inv;
fsum = fpair + fpair1;
// derivatives of the product of rij and ni, the result is a vector
dprodnorm1[0] =
dnormdri[i][0][0] * delx + dnormdri[i][1][0] * dely + dnormdri[i][2][0] * delz;
dprodnorm1[1] =
dnormdri[i][0][1] * delx + dnormdri[i][1][1] * dely + dnormdri[i][2][1] * delz;
dprodnorm1[2] =
dnormdri[i][0][2] * delx + dnormdri[i][1][2] * dely + dnormdri[i][2][2] * delz;
fp1[0] = prodnorm1 * normal[i][0] * fpair1;
fp1[1] = prodnorm1 * normal[i][1] * fpair1;
fp1[2] = prodnorm1 * normal[i][2] * fpair1;
fprod1[0] = prodnorm1 * dprodnorm1[0] * fpair1;
fprod1[1] = prodnorm1 * dprodnorm1[1] * fpair1;
fprod1[2] = prodnorm1 * dprodnorm1[2] * fpair1;
fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx / r;
fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely / r;
fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz / r;
f[i][0] += fkcx - fprod1[0] * Tap;
f[i][1] += fkcy - fprod1[1] * Tap;
f[i][2] += fkcz - fprod1[2] * Tap;
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
// calculate the forces acted on the neighbors of atom i from atom j
ILP_neighs_i = ILP_firstneigh[i];
for (kk = 0; kk < ILP_numneigh[i]; kk++) {
k = ILP_neighs_i[kk];
if (k == i) continue;
// derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
dprodnorm1[0] = dnormal[i][0][kk][0] * delx + dnormal[i][1][kk][0] * dely +
dnormal[i][2][kk][0] * delz;
dprodnorm1[1] = dnormal[i][0][kk][1] * delx + dnormal[i][1][kk][1] * dely +
dnormal[i][2][kk][1] * delz;
dprodnorm1[2] = dnormal[i][0][kk][2] * delx + dnormal[i][1][kk][2] * dely +
dnormal[i][2][kk][2] * delz;
fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap;
fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap;
fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap;
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
delki[0] = x[k][0] - x[i][0];
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag)
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz);
}
}
}
}
/* ----------------------------------------------------------------------
create ILP neighbor list from main neighbor list to calculate normals
------------------------------------------------------------------------- */
void PairILPTMD::ILP_neigh()
{
int i, j, l, ii, jj, ll, n, allnum, jnum, itype, jtype, ltype, imol, jmol, count;
double xtmp, ytmp, ztmp, delx, dely, delz, deljx, deljy, deljz, rsq, rsqlj;
int *ilist, *jlist, *numneigh, **firstneigh;
int *neighsort;
int neighptr[10], check[10];
double **x = atom->x;
int *type = atom->type;
if (atom->nmax > maxlocal) {
maxlocal = atom->nmax;
memory->destroy(ILP_numneigh);
memory->sfree(ILP_firstneigh);
memory->create(ILP_numneigh, maxlocal, "ILPTMD:numneigh");
ILP_firstneigh = (int **) memory->smalloc(maxlocal * sizeof(int *), "ILPTMD:firstneigh");
}
allnum = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// store all ILP neighs of owned and ghost atoms
// scan full neighbor list of I
ipage->reset();
for (ii = 0; ii < allnum; ii++) {
i = ilist[ii];
//initialize varibles
n = 0;
neighsort = ipage->vget();
for (ll = 0; ll < 10; ll++) {
neighptr[ll] = -1;
check[ll] = -1;
}
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = map[type[i]];
imol = atom->molecule[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
jmol = atom->molecule[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// check if the atom i is TMD, i.e., Mo/S/W/Se
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
strcmp(elements[itype], "Te") == 0) {
if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol && type[i] == type[j]) {
neighptr[n++] = j;
}
} else { // atom i is C, B, N or H.
if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol) { neighptr[n++] = j; }
}
} // loop over jj
// if atom i is Mo/W/S/Se/Te, then sorting the orders of neighbors
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 ||
strcmp(elements[itype], "Te") == 0) {
// initialize neighsort
for (ll = 0; ll < n; ll++) {
neighsort[ll] = neighptr[ll];
check[ll] = neighptr[ll];
}
// select the first neighbor of atomi
if (n == Nnei) {
neighsort[0] = neighptr[0];
check[0] = -1;
} else if (n < Nnei && n > 0) {
for (jj = 0; jj < n; jj++) { //identify the first neighbor
j = neighptr[jj];
jtype = map[type[j]];
count = 0;
for (ll = 0; ll < n; ll++) {
l = neighptr[ll];
ltype = map[type[l]];
if (l == j) continue;
deljx = x[l][0] - x[j][0];
deljy = x[l][1] - x[j][1];
deljz = x[l][2] - x[j][2];
rsqlj = deljx * deljx + deljy * deljy + deljz * deljz;
if (rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) { count++; }
}
if (count == 1) {
neighsort[0] = neighptr[jj];
check[jj] = -1;
break;
}
} // end of idenfying the first neighbor
} else if (n > Nnei) {
fprintf(screen, "Molecule ID = %d\n", imol);
fprintf(screen, "Atom Type = %d\n", type[i]);
fprintf(screen, "Neinum = %d\n", n);
error->one(FLERR,
"There are too many neighbors for TMD atoms, please check your configuration");
}
// sort the order of neighbors of atomi
for (jj = 0; jj < n; jj++) {
j = neighsort[jj];
jtype = map[type[j]];
ll = 0;
while (ll < n) {
l = neighptr[ll];
if (check[ll] == -1) {
ll++;
continue;
}
ltype = map[type[l]];
deljx = x[l][0] - x[j][0];
deljy = x[l][1] - x[j][1];
deljz = x[l][2] - x[j][2];
rsqlj = deljx * deljx + deljy * deljy + deljz * deljz;
if (rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) {
neighsort[jj + 1] = l;
check[ll] = -1;
break;
}
ll++;
}
} // end of sorting the order of neighbors
} else { // for B/N/C/H atoms
if (n > 3)
error->one(
FLERR,
"There are too many neighbors for B/N/C/H atoms, please check your configuration");
for (ll = 0; ll < n; ll++) { neighsort[ll] = neighptr[ll]; }
}
ILP_firstneigh[i] = neighsort;
ILP_numneigh[i] = n;
ipage->vgot(n);
if (ipage->status()) error->one(FLERR, "Neighbor list overflow, boost neigh_modify one");
}
}
/* ----------------------------------------------------------------------
Calculate the normals for each atom
------------------------------------------------------------------------- */
void PairILPTMD::calc_normal()
{
int i, j, ii, jj, inum, jnum;
int cont, id, ip, m, k, itype;
double nn, xtp, ytp, ztp, delx, dely, delz, nn2;
int *ilist, *jlist;
double Nave[3], dni[3], dpvdri[3][3];
double **x = atom->x;
int *type = atom->type;
memory->destroy(dnn);
memory->destroy(vect);
memory->destroy(pvet);
memory->destroy(dpvet1);
memory->destroy(dpvet2);
memory->destroy(dNave);
memory->create(dnn, Nnei, 3, "ILPTMD:dnn");
memory->create(vect, Nnei, 3, "ILPTMD:vect");
memory->create(pvet, Nnei, 3, "ILPTMD:pvet");
memory->create(dpvet1, Nnei, 3, 3, "ILPTMD:dpvet1");
memory->create(dpvet2, Nnei, 3, 3, "ILPTMD:dpvet2");
memory->create(dNave, 3, Nnei, 3, "ILPTMD:dNave");
// grow normal array if necessary
if (atom->nmax > nmax) {
memory->destroy(normal);
memory->destroy(dnormal);
memory->destroy(dnormdri);
nmax = atom->nmax;
memory->create(normal, nmax, 3, "ILPTMD:normal");
memory->create(dnormdri, nmax, 3, 3, "ILPTMD:dnormdri");
memory->create(dnormal, nmax, 3, Nnei, 3, "ILPTMD:dnormal");
}
inum = list->inum;
ilist = list->ilist;
//Calculate normals
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = map[type[i]];
xtp = x[i][0];
ytp = x[i][1];
ztp = x[i][2];
// Initialize the arrays
for (id = 0; id < 3; id++) {
Nave[id] = 0.0;
dni[id] = 0.0;
normal[i][id] = 0.0;
for (ip = 0; ip < 3; ip++) {
dpvdri[ip][id] = 0.0;
dnormdri[i][ip][id] = 0.0;
for (m = 0; m < Nnei; m++) {
dnn[m][id] = 0.0;
vect[m][id] = 0.0;
pvet[m][id] = 0.0;
dpvet1[m][ip][id] = 0.0;
dpvet2[m][ip][id] = 0.0;
dNave[id][m][ip] = 0.0;
dnormal[i][id][m][ip] = 0.0;
}
}
}
cont = 0;
jlist = ILP_firstneigh[i];
jnum = ILP_numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = x[j][0] - xtp;
dely = x[j][1] - ytp;
delz = x[j][2] - ztp;
vect[cont][0] = delx;
vect[cont][1] = dely;
vect[cont][2] = delz;
cont++;
}
//############################ For the dangling atoms ############################
if (cont <= 1) {
normal[i][0] = 0.0;
normal[i][1] = 0.0;
normal[i][2] = 1.0;
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = 0.0;
for (m = 0; m < Nnei; m++) { dnormal[i][id][m][ip] = 0.0; }
}
}
}
//############################ For the edge atoms of TMD ################################
else if (cont > 1 && cont < Nnei) {
if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 ||
strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0) {
// derivatives of Ni[l] respect to the cont neighbors
for (k = 0; k < cont - 1; k++) {
for (ip = 0; ip < 3; ip++) {
pvet[k][ip] = vect[k][modulo(ip + 1, 3)] * vect[k + 1][modulo(ip + 2, 3)] -
vect[k][modulo(ip + 2, 3)] * vect[k + 1][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[k][0][0] = 0.0;
dpvet1[k][0][1] = vect[modulo(k + 1, Nnei)][2];
dpvet1[k][0][2] = -vect[modulo(k + 1, Nnei)][1];
// dNik,y/drl
dpvet1[k][1][0] = -vect[modulo(k + 1, Nnei)][2];
dpvet1[k][1][1] = 0.0;
dpvet1[k][1][2] = vect[modulo(k + 1, Nnei)][0];
// dNik,z/drl
dpvet1[k][2][0] = vect[modulo(k + 1, Nnei)][1];
dpvet1[k][2][1] = -vect[modulo(k + 1, Nnei)][0];
dpvet1[k][2][2] = 0.0;
// dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNik,x/drl+1
dpvet2[k][0][0] = 0.0;
dpvet2[k][0][1] = -vect[modulo(k, Nnei)][2];
dpvet2[k][0][2] = vect[modulo(k, Nnei)][1];
// dNik,y/drl+1
dpvet2[k][1][0] = vect[modulo(k, Nnei)][2];
dpvet2[k][1][1] = 0.0;
dpvet2[k][1][2] = -vect[modulo(k, Nnei)][0];
// dNik,z/drl+1
dpvet2[k][2][0] = -vect[modulo(k, Nnei)][1];
dpvet2[k][2][1] = vect[modulo(k, Nnei)][0];
dpvet2[k][2][2] = 0.0;
}
// average the normal vectors by using the Nnei neighboring planes
for (ip = 0; ip < 3; ip++) {
Nave[ip] = 0.0;
for (k = 0; k < cont - 1; k++) { Nave[ip] += pvet[k][ip]; }
Nave[ip] /= (cont - 1);
}
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// derivatives of non-normalized normal vector, dNave:3xcontx3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (m = 0; m < cont; m++) {
if (m == 0) {
dNave[id][m][ip] = dpvet1[m][id][ip] / (cont - 1);
} else if (m == cont - 1) {
dNave[id][m][ip] = dpvet2[m - 1][id][ip] / (cont - 1);
} else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m - 1][id][ip]) / (cont - 1);
}
}
}
}
// derivatives of nn, dnn:contx3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
// Calculte dNave/dri, defined as dpvdri
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dpvdri[id][ip] = 0.0;
for (k = 0; k < cont; k++) { dpvdri[id][ip] -= dNave[id][k][ip]; }
}
}
// derivatives of nn, dnn:3x1 vector
dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn;
dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn;
dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn;
// derivatives of unit vector ni respect to ri, the result is 3x3 matrix
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2;
}
}
} // for TMD
//############################ For the edge & bulk atoms of GrhBN ################################
else {
if (cont == 2) {
for (ip = 0; ip < 3; ip++) {
pvet[0][ip] = vect[0][modulo(ip + 1, 3)] * vect[1][modulo(ip + 2, 3)] -
vect[0][modulo(ip + 2, 3)] * vect[1][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[0][0][0] = 0.0;
dpvet1[0][0][1] = vect[1][2];
dpvet1[0][0][2] = -vect[1][1];
// dNi0,y/drl
dpvet1[0][1][0] = -vect[1][2];
dpvet1[0][1][1] = 0.0;
dpvet1[0][1][2] = vect[1][0];
// dNi0,z/drl
dpvet1[0][2][0] = vect[1][1];
dpvet1[0][2][1] = -vect[1][0];
dpvet1[0][2][2] = 0.0;
// dpvet2[0][l][ip]: the derivatve of the 0 (=0,...cont-1)th Ni0 respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNi0,x/drl+1
dpvet2[0][0][0] = 0.0;
dpvet2[0][0][1] = -vect[0][2];
dpvet2[0][0][2] = vect[0][1];
// dNi0,y/drl+1
dpvet2[0][1][0] = vect[0][2];
dpvet2[0][1][1] = 0.0;
dpvet2[0][1][2] = -vect[0][0];
// dNi0,z/drl+1
dpvet2[0][2][0] = -vect[0][1];
dpvet2[0][2][1] = vect[0][0];
dpvet2[0][2][2] = 0.0;
for (ip = 0; ip < 3; ip++) { Nave[ip] += pvet[0][ip]; }
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// derivatives of non-normalized normal vector, dNave:3xcontx3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (m = 0; m < cont; m++) {
if (m == 0) {
dNave[id][m][ip] = dpvet1[m][id][ip] / (cont - 1);
} else if (m == cont - 1) {
dNave[id][m][ip] = dpvet2[m - 1][id][ip] / (cont - 1);
} else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m - 1][id][ip]) / (cont - 1);
}
}
}
}
// derivatives of nn, dnn:contx3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1
for (m = 0; m < cont; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
// Calculte dNave/dri, defined as dpvdri
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dpvdri[id][ip] = 0.0;
for (k = 0; k < cont; k++) { dpvdri[id][ip] -= dNave[id][k][ip]; }
}
}
// derivatives of nn, dnn:3x1 vector
dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn;
dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn;
dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn;
// derivatives of unit vector ni respect to ri, the result is 3x3 matrix
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2;
}
}
} // end of cont == 2
else if (cont == 3) {
// derivatives of Ni[l] respect to the 3 neighbors
for (k = 0; k < 3; k++) {
for (ip = 0; ip < 3; ip++) {
pvet[k][ip] = vect[modulo(k, 3)][modulo(ip + 1, 3)] *
vect[modulo(k + 1, 3)][modulo(ip + 2, 3)] -
vect[modulo(k, 3)][modulo(ip + 2, 3)] * vect[modulo(k + 1, 3)][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[k][0][0] = 0.0;
dpvet1[k][0][1] = vect[modulo(k + 1, 3)][2];
dpvet1[k][0][2] = -vect[modulo(k + 1, 3)][1];
// dNik,y/drl
dpvet1[k][1][0] = -vect[modulo(k + 1, 3)][2];
dpvet1[k][1][1] = 0.0;
dpvet1[k][1][2] = vect[modulo(k + 1, 3)][0];
// dNik,z/drl
dpvet1[k][2][0] = vect[modulo(k + 1, 3)][1];
dpvet1[k][2][1] = -vect[modulo(k + 1, 3)][0];
dpvet1[k][2][2] = 0.0;
// dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNik,x/drl+1
dpvet2[k][0][0] = 0.0;
dpvet2[k][0][1] = -vect[modulo(k, 3)][2];
dpvet2[k][0][2] = vect[modulo(k, 3)][1];
// dNik,y/drl+1
dpvet2[k][1][0] = vect[modulo(k, 3)][2];
dpvet2[k][1][1] = 0.0;
dpvet2[k][1][2] = -vect[modulo(k, 3)][0];
// dNik,z/drl+1
dpvet2[k][2][0] = -vect[modulo(k, 3)][1];
dpvet2[k][2][1] = vect[modulo(k, 3)][0];
dpvet2[k][2][2] = 0.0;
}
// average the normal vectors by using the 3 neighboring planes
for (ip = 0; ip < 3; ip++) {
Nave[ip] = 0.0;
for (k = 0; k < 3; k++) { Nave[ip] += pvet[k][ip]; }
Nave[ip] /= 3;
}
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// for the central atoms, dnormdri is always zero
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) { dnormdri[i][id][ip] = 0.0; }
}
// derivatives of non-normalized normal vector, dNave:3x3x3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (
m = 0; m < 3;
m++) { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] =
(dpvet1[modulo(m, 3)][id][ip] + dpvet2[modulo(m - 1, 3)][id][ip]) / 3;
}
}
}
// derivatives of nn, dnn:3x3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...3-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < 3; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] +
Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,3-1
for (m = 0; m < 3; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
} // end of cont == 3
else
error->one(FLERR,
"There are too many neighbors for calculating normals of B/N/C/H atoms");
} // for B/N/C/H
} // end of if(cont<Nnei)
//############################ For the bulk atoms ################################
else if (cont == Nnei) {
// derivatives of Ni[l] respect to the Nnei neighbors
for (k = 0; k < Nnei; k++) {
for (ip = 0; ip < 3; ip++) {
pvet[k][ip] = vect[modulo(k, Nnei)][modulo(ip + 1, 3)] *
vect[modulo(k + 1, Nnei)][modulo(ip + 2, 3)] -
vect[modulo(k, Nnei)][modulo(ip + 2, 3)] *
vect[modulo(k + 1, Nnei)][modulo(ip + 1, 3)];
}
// dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l
// derivatives respect to atom l
// dNik,x/drl
dpvet1[k][0][0] = 0.0;
dpvet1[k][0][1] = vect[modulo(k + 1, Nnei)][2];
dpvet1[k][0][2] = -vect[modulo(k + 1, Nnei)][1];
// dNik,y/drl
dpvet1[k][1][0] = -vect[modulo(k + 1, Nnei)][2];
dpvet1[k][1][1] = 0.0;
dpvet1[k][1][2] = vect[modulo(k + 1, Nnei)][0];
// dNik,z/drl
dpvet1[k][2][0] = vect[modulo(k + 1, Nnei)][1];
dpvet1[k][2][1] = -vect[modulo(k + 1, Nnei)][0];
dpvet1[k][2][2] = 0.0;
// dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1
// derivatives respect to atom l+1
// dNik,x/drl+1
dpvet2[k][0][0] = 0.0;
dpvet2[k][0][1] = -vect[modulo(k, Nnei)][2];
dpvet2[k][0][2] = vect[modulo(k, Nnei)][1];
// dNik,y/drl+1
dpvet2[k][1][0] = vect[modulo(k, Nnei)][2];
dpvet2[k][1][1] = 0.0;
dpvet2[k][1][2] = -vect[modulo(k, Nnei)][0];
// dNik,z/drl+1
dpvet2[k][2][0] = -vect[modulo(k, Nnei)][1];
dpvet2[k][2][1] = vect[modulo(k, Nnei)][0];
dpvet2[k][2][2] = 0.0;
}
// average the normal vectors by using the Nnei neighboring planes
for (ip = 0; ip < 3; ip++) {
Nave[ip] = 0.0;
for (k = 0; k < Nnei; k++) { Nave[ip] += pvet[k][ip]; }
Nave[ip] /= Nnei;
}
// the magnitude of the normal vector
nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2];
nn = sqrt(nn2);
if (nn == 0.0) error->one(FLERR, "The magnitude of the normal vector is zero");
// the unit normal vector
normal[i][0] = Nave[0] / nn;
normal[i][1] = Nave[1] / nn;
normal[i][2] = Nave[2] / nn;
// for the central atoms, dnormdri is always zero
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) { dnormdri[i][id][ip] = 0.0; }
}
// derivatives of non-normalized normal vector, dNave:3xNneix3 array
// dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
for (
m = 0; m < Nnei;
m++) { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m
dNave[id][m][ip] =
(dpvet1[modulo(m, Nnei)][id][ip] + dpvet2[modulo(m - 1, Nnei)][id][ip]) / Nnei;
}
}
}
// derivatives of nn, dnn:Nneix3 vector
// dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2
// r[m][id]: the id's component of atom m
for (m = 0; m < Nnei; m++) {
for (id = 0; id < 3; id++) {
dnn[m][id] =
(Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] + Nave[2] * dNave[2][m][id]) /
nn;
}
}
// dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2.
// for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1
for (m = 0; m < Nnei; m++) {
for (id = 0; id < 3; id++) {
for (ip = 0; ip < 3; ip++) {
dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2;
}
}
}
} else {
error->one(FLERR, "There are too many neighbors for calculating normals of TMD atoms");
} // end of four cases of cont
} // end of i loop
}

View File

@ -0,0 +1,70 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(ilp/tmd,PairILPTMD);
// clang-format on
#else
#ifndef LMP_PAIR_ILP_TMD_H
#define LMP_PAIR_ILP_TMD_H
#include "pair_ilp_graphene_hbn.h"
namespace LAMMPS_NS {
class PairILPTMD : public PairILPGrapheneHBN {
public:
PairILPTMD(class LAMMPS *);
protected:
void settings(int, char **) override;
void ILP_neigh() override;
void calc_normal() override;
void calc_FRep(int, int) override;
/**************************************************************/
/* modulo operation with cycling around range */
inline int modulo(int k, int range)
{
if (k < 0) k += range;
return k % range;
}
/**************************************************************/
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/

View File

@ -28,8 +28,8 @@
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include <cmath>
@ -160,15 +160,15 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
if (evflag) {
ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
if (vflag_either) {
double fi[3],fj[3];
double fi[3], fj[3];
fi[0] = delx * fpair1;
fi[1] = dely * fpair1;
fi[2] = 0;
fj[0] = -delx * fpair1;
fj[1] = -dely * fpair1;
fj[2] = 0;
v_tally2_newton(i,fi,x[i]);
v_tally2_newton(j,fj,x[j]);
v_tally2_newton(i, fi, x[i]);
v_tally2_newton(j, fj, x[j]);
}
}
}
@ -246,9 +246,9 @@ void PairKolmogorovCrespiZ::coeff(int narg, char **arg)
void PairKolmogorovCrespiZ::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style kolmogorov/crespi/z requires newton pair on");
error->all(FLERR, "Pair style kolmogorov/crespi/z requires newton pair on");
neighbor->request(this,instance_me);
neighbor->request(this, instance_me);
}
/* ----------------------------------------------------------------------

View File

@ -0,0 +1,245 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
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: Wengen Ouyang (Wuhan University)
e-mail: w.g.ouyang at gmail dot com
This is a full version of the potential described in
[Ouyang et al, J. Chem. Theory Comput. 17, 7215-7223 (2021)]
------------------------------------------------------------------------- */
#include "pair_saip_metal.h"
#include "atom.h"
#include "citeme.h"
#include "error.h"
#include "force.h"
#include "interlayer_taper.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace InterLayer;
#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
static const char cite_saip[] =
"saip/metal potential doi.org/10.1021/acs.jctc.1c00622\n"
"@Article{Ouyang2021\n"
" author = {W. Ouyang, O. Hod, and R. Guerra},\n"
" title = {Registry-Dependent Potential for Interfaces of Gold with Graphitic Systems},\n"
" journal = {J. Chem. Theory Comput.},\n"
" volume = 17,\n"
" pages = {7215-7223}\n"
" year = 2021,\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairSAIPMETAL::PairSAIPMETAL(LAMMPS *lmp) : PairILPGrapheneHBN(lmp)
{
variant = SAIP_METAL;
single_enable = 0;
if (lmp->citeme) lmp->citeme->add(cite_saip);
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSAIPMETAL::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR, "Illegal pair_style command");
if (!utils::strmatch(force->pair_style, "^hybrid/overlay"))
error->all(FLERR, "Pair style saip/metal must be used as sub-style with hybrid/overlay");
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp);
}
/* ----------------------------------------------------------------------
Repulsive forces and energy
------------------------------------------------------------------------- */
void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */)
{
int i, j, ii, jj, inum, jnum, itype, jtype, k, kk;
double prodnorm1, fkcx, fkcy, fkcz, filp;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1;
double rsq, r, Rcut, rhosq1, exp0, exp1, Tap, dTap, Vilp;
double frho1, Erep, fsum, rdsq1;
int *ilist, *jlist, *numneigh, **firstneigh;
int *ILP_neighs_i;
evdwl = 0.0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double dprodnorm1[3] = {0.0, 0.0, 0.0};
double fp1[3] = {0.0, 0.0, 0.0};
double fprod1[3] = {0.0, 0.0, 0.0};
double delki[3] = {0.0, 0.0, 0.0};
double fk[3] = {0.0, 0.0, 0.0};
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
//calculate exp(-lambda*(r-z0))*[epsilon/2 + f(rho_ij)]
// loop over neighbors of owned atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// only include the interaction between different layers
if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) {
int iparam_ij = elem2param[map[itype]][map[jtype]];
Param &p = params[iparam_ij];
r = sqrt(rsq);
// turn on/off taper function
if (tap_flag) {
Rcut = sqrt(cutsq[itype][jtype]);
Tap = calc_Tap(r, Rcut);
dTap = calc_dTap(r, Rcut);
} else {
Tap = 1.0;
dTap = 0.0;
}
// for atoms in bulk materials
if (strcmp(elements[map[itype]], "C") != 0 && strcmp(elements[map[itype]], "H") != 0 &&
strcmp(elements[map[itype]], "B") != 0 && strcmp(elements[map[itype]], "N") != 0) {
// Set ni along rij
exp0 = exp(-p.lambda * (r - p.z0));
frho1 = 1.0 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda * exp0 / r * Erep;
filp = fpair * Tap - Vilp * dTap / r;
fkcx = delx * filp;
fkcy = dely * filp;
fkcz = delz * filp;
f[i][0] += fkcx;
f[i][1] += fkcy;
f[i][2] += fkcz;
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, filp, delx, dely, delz);
} else { // for atoms in 2D materials
// Calculate the transverse distance
prodnorm1 = normal[i][0] * delx + normal[i][1] * dely + normal[i][2] * delz;
rhosq1 = rsq - prodnorm1 * prodnorm1; // rho_ij
rdsq1 = rhosq1 * p.delta2inv; // (rho_ij/delta)^2
// store exponents
exp0 = exp(-p.lambda * (r - p.z0));
exp1 = exp(-rdsq1);
frho1 = exp1 * p.C;
Erep = 0.5 * p.epsilon + frho1;
Vilp = exp0 * Erep;
// derivatives
fpair = p.lambda * exp0 / r * Erep;
fpair1 = 2.0 * exp0 * frho1 * p.delta2inv;
fsum = fpair + fpair1;
// derivatives of the product of rij and ni, the result is a vector
dprodnorm1[0] =
dnormdri[0][0][i] * delx + dnormdri[1][0][i] * dely + dnormdri[2][0][i] * delz;
dprodnorm1[1] =
dnormdri[0][1][i] * delx + dnormdri[1][1][i] * dely + dnormdri[2][1][i] * delz;
dprodnorm1[2] =
dnormdri[0][2][i] * delx + dnormdri[1][2][i] * dely + dnormdri[2][2][i] * delz;
fp1[0] = prodnorm1 * normal[i][0] * fpair1;
fp1[1] = prodnorm1 * normal[i][1] * fpair1;
fp1[2] = prodnorm1 * normal[i][2] * fpair1;
fprod1[0] = prodnorm1 * dprodnorm1[0] * fpair1;
fprod1[1] = prodnorm1 * dprodnorm1[1] * fpair1;
fprod1[2] = prodnorm1 * dprodnorm1[2] * fpair1;
fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx / r;
fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely / r;
fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz / r;
f[i][0] += fkcx - fprod1[0] * Tap;
f[i][1] += fkcy - fprod1[1] * Tap;
f[i][2] += fkcz - fprod1[2] * Tap;
f[j][0] -= fkcx;
f[j][1] -= fkcy;
f[j][2] -= fkcz;
// calculate the forces acted on the neighbors of atom i from atom j
ILP_neighs_i = ILP_firstneigh[i];
for (kk = 0; kk < ILP_numneigh[i]; kk++) {
k = ILP_neighs_i[kk];
if (k == i) continue;
// derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i
dprodnorm1[0] = dnormal[0][0][kk][i] * delx + dnormal[1][0][kk][i] * dely +
dnormal[2][0][kk][i] * delz;
dprodnorm1[1] = dnormal[0][1][kk][i] * delx + dnormal[1][1][kk][i] * dely +
dnormal[2][1][kk][i] * delz;
dprodnorm1[2] = dnormal[0][2][kk][i] * delx + dnormal[1][2][kk][i] * dely +
dnormal[2][2][kk][i] * delz;
fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap;
fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap;
fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap;
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
delki[0] = x[k][0] - x[i][0];
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}
if (eflag) pvector[1] += evdwl = Tap * Vilp;
if (evflag)
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz);
} // end of speration atoms for bulk materials
} // end of rsq < cutoff
} // loop over jj
} // loop over ii
}

View File

@ -0,0 +1,58 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(saip/metal,PairSAIPMETAL);
// clang-format on
#else
#ifndef LMP_PAIR_SAIP_METAL_H
#define LMP_PAIR_SAIP_METAL_H
#include "pair_ilp_graphene_hbn.h"
namespace LAMMPS_NS {
class PairSAIPMETAL : public PairILPGrapheneHBN {
public:
PairSAIPMETAL(class LAMMPS *);
protected:
void settings(int, char **) override;
void calc_FRep(int, int) override;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/

View File

@ -207,6 +207,11 @@ KokkosLMP::KokkosLMP(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
error->all(FLERR,"Kokkos has been compiled with GPU-enabled backend but no GPUs are requested");
#endif
#if !defined(KOKKOS_ENABLE_OPENMP) && !defined(KOKKOS_ENABLE_THREADS)
if (nthreads > 1)
error->all(FLERR,"Multiple CPU threads are requested but Kokkos has not been compiled using a threading-enabled backend");
#endif
#ifndef KOKKOS_ENABLE_SERIAL
if (nthreads == 1 && me == 0)
error->warning(FLERR,"When using a single thread, the Kokkos Serial backend "

View File

@ -104,6 +104,15 @@ E: Kokkos has been compiled with GPU-enabled backend but no GPUs are requested
One or more GPUs must be used when Kokkos is compiled for CUDA/HIP/SYCL/OpenMPTarget.
E: Multiple CPU threads are requested but Kokkos has not been compiled using a threading-enabled backend
Must use the Kokkos OpenMP or Threads backend for multiple threads.
W: When using a single thread, the Kokkos Serial backend (i.e. Makefile.kokkos_mpi_only)
gives better performance than the OpenMP backend
Self-expanatory.
W: Kokkos package already initalized, cannot reinitialize with different parameters
Self-explanatory.

View File

@ -579,6 +579,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
FixBondReact::~FixBondReact()
{
for (int i = 0; i < narrhenius; i++) {
delete rrhandom[i];
}
delete [] rrhandom;
for (int i = 0; i < nreacts; i++) {
delete random[i];
}

View File

@ -68,14 +68,14 @@ void PairHybridScaled::compute(int eflag, int vflag)
const int nvars = scalevars.size();
if (nvars > 0) {
double *vals = new double[nvars];
for (i = 0; i < nvars; ++i) {
j = input->variable->find(scalevars[i].c_str());
if (j < 0)
error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[i]);
vals[i] = input->variable->compute_equal(j);
for (int k = 0; k < nvars; ++k) {
int m = input->variable->find(scalevars[k].c_str());
if (m < 0)
error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]);
vals[k] = input->variable->compute_equal(m);
}
for (i = 0; i < nstyles; ++i) {
if (scaleidx[i] >= 0) scaleval[i] = vals[scaleidx[i]];
for (int k = 0; k < nstyles; ++k) {
if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]];
}
delete[] vals;
}
@ -386,14 +386,14 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq,
const int nvars = scalevars.size();
if (nvars > 0) {
double *vals = new double[nvars];
for (i = 0; i < nvars; ++i) {
j = input->variable->find(scalevars[i].c_str());
if (j < 0)
error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[i]);
vals[i] = input->variable->compute_equal(j);
for (int k = 0; k < nvars; ++k) {
int m = input->variable->find(scalevars[k].c_str());
if (m < 0)
error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]);
vals[k] = input->variable->compute_equal(m);
}
for (i = 0; i < nstyles; ++i) {
if (scaleidx[i] >= 0) scaleval[i] = vals[scaleidx[i]];
for (int k = 0; k < nstyles; ++k) {
if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]];
}
delete[] vals;
}
@ -593,7 +593,7 @@ void PairHybridScaled::init_svector()
void PairHybridScaled::copy_svector(int itype, int jtype)
{
int n = 0;
Pair *this_style;
Pair *this_style = nullptr;
// fill svector array.
// copy data from active styles and use 0.0 for inactive ones

View File

@ -1 +1 @@
#define LAMMPS_VERSION "7 Jan 2022"
#define LAMMPS_VERSION "17 Feb 2022"

View File

@ -236,7 +236,7 @@ TEST_F(LAMMPS_omp, InitMembers)
}
}
// test fixture for Kokkos tests
// test fixture for Kokkos using 2 threads if threading is available
class LAMMPS_kokkos : public ::testing::Test {
protected:
LAMMPS *lmp;
@ -256,18 +256,15 @@ protected:
void SetUp() override
{
const char *args[] = {"LAMMPS_test", "-log", "none", "-echo", "none", "-screen", "none",
"-k", "on", "t", "2", "-sf", "kk"};
char **argv = (char **)args;
int argc = sizeof(args) / sizeof(char *);
// only run this test fixture with kk suffix if KOKKOS package is installed
"-k", "on", "t", "1", "-sf", "kk"};
if (Info::has_accelerator_feature("KOKKOS", "api", "openmp")) args[10] = "2";
char **argv = (char **)args;
int argc = sizeof(args) / sizeof(char *);
if (LAMMPS::is_installed_pkg("KOKKOS")) {
::testing::internal::CaptureStdout();
lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
std::string output = ::testing::internal::GetCapturedStdout();
if (Info::has_accelerator_feature("KOKKOS", "api", "openmp"))
EXPECT_THAT(output, StartsWith("Kokkos::OpenMP::"));
lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
::testing::internal::GetCapturedStdout();
} else
GTEST_SKIP();
}
@ -307,6 +304,10 @@ TEST_F(LAMMPS_kokkos, InitMembers)
EXPECT_EQ(lmp->num_package, 0);
EXPECT_EQ(lmp->clientserver, 0);
if (Info::has_accelerator_feature("KOKKOS", "api", "openmp"))
EXPECT_EQ(lmp->comm->nthreads, 2);
else
EXPECT_EQ(lmp->comm->nthreads, 1);
EXPECT_NE(lmp->kokkos, nullptr);
EXPECT_NE(lmp->atomKK, nullptr);
EXPECT_NE(lmp->memoryKK, nullptr);

View File

@ -1,7 +1,7 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:09:03 2021
epsilon: 5e-13
epsilon: 1e-12
prerequisites: ! |
pair kim
pre_commands: ! |

View File

@ -1,7 +1,7 @@
---
lammps_version: 30 Jul 2021
date_generated: Wed Aug 25 07:37:07 2021
epsilon: 2e-12
epsilon: 2e-9
skip_tests: single
prerequisites: ! |
pair lebedeva/z

View File

@ -1,7 +1,7 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:08:41 2021
epsilon: 2e-08
epsilon: 5e-08
prerequisites: ! |
atom full
pair buck/long/coul/long

View File

@ -1,7 +1,7 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:08:53 2021
epsilon: 2e-08
epsilon: 5e-08
prerequisites: ! |
atom full
pair lj/long/coul/long

View File

@ -1,6 +1,7 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:08:56 2021
tags: unstable
epsilon: 1.5e-12
prerequisites: ! |
atom full