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

start document describing the implementation of a pair style

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-02-27 01:58:59 -05:00
parent c8f3f4b173
commit 93b42f14e9
4 changed files with 560 additions and 245 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Build_development.rst
View File

@ -523,6 +523,8 @@ The following options are available.
These should help to make source and documentation files conforming
to some the coding style preferences of the LAMMPS developers.
.. _clang-format:
Clang-format support
--------------------

249
doc/src/Developer_write.rst
View File

@ -6,250 +6,9 @@ be extended by writing new classes that derive from existing
parent classes in LAMMPS. Here, some specific coding
details are provided for writing code for LAMMPS.
Writing a new fix style
^^^^^^^^^^^^^^^^^^^^^^^
Writing fixes is a flexible way of extending LAMMPS. Users can
implement many things using fixes:
- changing particles attributes (positions, velocities, forces, etc.). Examples: FixNVE, FixFreeze.
- reading/writing data. Example: FixRestart.
- adding or modifying properties due to geometry. Example: FixWall.
- interacting with other subsystems or external code: Examples: FixTTM, FixExternal, FixLATTE
- saving information for analysis or future use (previous positions,
for instance). Examples: Fix AveTime, FixStoreState.
All fixes are derived from the Fix base class and must have a
constructor with the signature: ``FixPrintVel(class LAMMPS *, int, char **)``.
Every fix must be registered in LAMMPS by writing the following lines
of code in the header before include guards:
.. code-block:: c++
#ifdef FIX_CLASS
// clang-format off
FixStyle(print/vel,FixPrintVel);
// clang-format on
#else
/* the definition of the FixPrintVel class comes here */
...
#endif
Where ``print/vel`` is the style name of your fix in the input script and
``FixPrintVel`` is the name of the class. The header file would be called
``fix_print_vel.h`` and the implementation file ``fix_print_vel.cpp``.
These conventions allow LAMMPS to automatically integrate it into the
executable when compiling and associate your new fix class with the designated
keyword when it parses the input script.
Let's write a simple fix which will print the average velocity at the end
of each timestep. First of all, implement a constructor:
.. code-block:: c++
FixPrintVel::FixPrintVel(LAMMPS *lmp, int narg, char **arg)
: Fix(lmp, narg, arg)
{
if (narg < 4)
error->all(FLERR,"Illegal fix print/vel command");
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery <= 0)
error->all(FLERR,"Illegal fix print/vel command");
}
In the constructor you should parse your fix arguments which are
specified in the script. All fixes have pretty much the same syntax:
``fix <fix-ID> <fix group> <fix name> <fix arguments ...>``. The
first 3 parameters are parsed by Fix base class constructor, while
``<fix arguments>`` should be parsed by you. In our case, we need to
specify how often we want to print an average velocity. For instance,
once in 50 timesteps: ``fix 1 print/vel 50``. There is a special variable
in the Fix class called ``nevery`` which specifies how often the method
``end_of_step()`` is called. Thus all we need to do is just set it up.
The next method we need to implement is ``setmask()``:
.. code-block:: c++
int FixPrintVel::setmask()
{
int mask = 0;
mask |= FixConst::END_OF_STEP;
return mask;
}
Here the user specifies which methods of your fix should be called
during execution. The constant ``END_OF_STEP`` corresponds to the
``end_of_step()`` method. The most important available methods that
are called during a timestep and the order in which they are called
are shown in the previous section.
.. code-block:: c++
void FixPrintVel::end_of_step()
{
// for add3, scale3
using namespace MathExtra;
double** v = atom->v;
int nlocal = atom->nlocal;
double localAvgVel[4]; // 4th element for particles count
memset(localAvgVel, 0, 4 * sizeof(double));
for (int particleInd = 0; particleInd < nlocal; ++particleInd) {
add3(localAvgVel, v[particleInd], localAvgVel);
}
localAvgVel[3] = nlocal;
double globalAvgVel[4];
memset(globalAvgVel, 0, 4 * sizeof(double));
MPI_Allreduce(localAvgVel, globalAvgVel, 4, MPI_DOUBLE, MPI_SUM, world);
scale3(1.0 / globalAvgVel[3], globalAvgVel);
if ((comm->me == 0) && screen) {
fmt::print(screen,"{}, {}, {}\n",
globalAvgVel[0], globalAvgVel[1], globalAvgVel[2]);
}
}
In the code above, we use MathExtra routines defined in
``math_extra.h``. There are bunch of math functions to work with
arrays of doubles as with math vectors. It is also important to note
that LAMMPS code should always assume to be run in parallel and that
atom data is thus distributed across the MPI ranks. Thus you can
only process data from local atoms directly and need to use MPI library
calls to combine or exchange data. For serial execution, LAMMPS
comes bundled with the MPI STUBS library that contains the MPI library
function calls in dummy versions that only work for a single MPI rank.
In this code we use an instance of Atom class. This object is stored
in the Pointers class (see ``pointers.h``) which is the base class of
the Fix base class. This object contains references to various class
instances (the original instances are created and held by the LAMMPS
class) with all global information about the simulation system.
Data from the Pointers class is available to all classes inherited from
it using protected inheritance. Hence when you write you own class,
which is going to use LAMMPS data, don't forget to inherit from Pointers
or pass an Pointer to it to all functions that need access. When writing
fixes we inherit from class Fix which is inherited from Pointers so
there is no need to inherit from it directly.
The code above computes average velocity for all particles in the
simulation. Yet you have one unused parameter in fix call from the
script: ``group_name``. This parameter specifies the group of atoms
used in the fix. So we should compute average for all particles in the
simulation only if ``group_name == "all"``, but it can be any group.
The group membership information of an atom is contained in the *mask*
property of and atom and the bit corresponding to a given group is
stored in the groupbit variable which is defined in Fix base class:
.. code-block:: c++
for (int i = 0; i < nlocal; ++i) {
if (atom->mask[i] & groupbit) {
// Do all job here
}
}
Class Atom encapsulates atoms positions, velocities, forces, etc. User
can access them using particle index. Note, that particle indexes are
usually changed every few timesteps because of neighbor list rebuilds
and spatial sorting (to improve cache efficiency).
Let us consider another Fix example: We want to have a fix which stores
atoms position from previous time step in your fix. The local atoms
indexes may not be valid on the next iteration. In order to handle
this situation there are several methods which should be implemented:
- ``double memory_usage()``: return how much memory the fix uses (optional)
- ``void grow_arrays(int)``: do reallocation of the per particle arrays in your fix
- ``void copy_arrays(int i, int j, int delflag)``: copy i-th per-particle
information to j-th. Used when atom sorting is performed. if delflag is set
and atom j owns a body, move the body information to atom i.
- ``void set_arrays(int i)``: sets i-th particle related information to zero
Note, that if your class implements these methods, it must call add calls of
add_callback and delete_callback to constructor and destructor. Since we want
to store positions of atoms from previous timestep, we need to add
``double** xold`` to the header file. Than add allocation code
to the constructor:
.. code-block:: c++
FixSavePos::FixSavePos(LAMMPS *lmp, int narg, char **arg), xold(nullptr)
{
//...
memory->create(xold, atom->nmax, 3, "FixSavePos:x");
atom->add_callback(0);
}
FixSavePos::~FixSavePos() {
atom->delete_callback(id, 0);
memory->destroy(xold);
}
Implement the aforementioned methods:
.. code-block:: c++
double FixSavePos::memory_usage()
{
int nmax = atom->nmax;
double bytes = 0.0;
bytes += nmax * 3 * sizeof(double);
return bytes;
}
void FixSavePos::grow_arrays(int nmax)
{
memory->grow(xold, nmax, 3, "FixSavePos:xold");
}
void FixSavePos::copy_arrays(int i, int j, int delflag)
{
memcpy(xold[j], xold[i], sizeof(double) * 3);
}
void FixSavePos::set_arrays(int i)
{
memset(xold[i], 0, sizeof(double) * 3);
}
int FixSavePos::pack_exchange(int i, double *buf)
{
int m = 0;
buf[m++] = xold[i][0];
buf[m++] = xold[i][1];
buf[m++] = xold[i][2];
return m;
}
int FixSavePos::unpack_exchange(int nlocal, double *buf)
{
int m = 0;
xold[nlocal][0] = buf[m++];
xold[nlocal][1] = buf[m++];
xold[nlocal][2] = buf[m++];
return m;
}
Now, a little bit about memory allocation. We use the Memory class which
is just a bunch of template functions for allocating 1D and 2D
arrays. So you need to add include ``memory.h`` to have access to them.
Finally, if you need to write/read some global information used in
your fix to the restart file, you might do it by setting flag
``restart_global = 1`` in the constructor and implementing methods void
``write_restart(FILE *fp)`` and ``void restart(char *buf)``.
If, in addition, you want to write the per-atom property to restart
files additional settings and functions are needed:
- a fix flag indicating this needs to be set ``restart_peratom = 1;``
- ``atom->add_callback()`` and ``atom->delete_callback()`` must be called
a second time with the final argument set to 1 instead of 0 (indicating
restart processing instead of per-atom data memory management).
- the functions ``void pack_restart(int i, double *buf)`` and
``void unpack_restart(int nlocal, int nth)`` need to be implemented
.. toctree::
:maxdepth: 1
Developer_write_pair
Developer_write_fix

246
doc/src/Developer_write_fix.rst Normal file
View File

@ -0,0 +1,246 @@
Writing a new fix style
^^^^^^^^^^^^^^^^^^^^^^^
Writing fixes is a flexible way of extending LAMMPS. Users can
implement many things using fixes:
- changing particles attributes (positions, velocities, forces, etc.). Examples: FixNVE, FixFreeze.
- reading/writing data. Example: FixRestart.
- adding or modifying properties due to geometry. Example: FixWall.
- interacting with other subsystems or external code: Examples: FixTTM, FixExternal, FixLATTE
- saving information for analysis or future use (previous positions,
for instance). Examples: Fix AveTime, FixStoreState.
All fixes are derived from the Fix base class and must have a
constructor with the signature: ``FixPrintVel(class LAMMPS *, int, char **)``.
Every fix must be registered in LAMMPS by writing the following lines
of code in the header before include guards:
.. code-block:: c++
#ifdef FIX_CLASS
// clang-format off
FixStyle(print/vel,FixPrintVel);
// clang-format on
#else
/* the definition of the FixPrintVel class comes here */
...
#endif
Where ``print/vel`` is the style name of your fix in the input script and
``FixPrintVel`` is the name of the class. The header file would be called
``fix_print_vel.h`` and the implementation file ``fix_print_vel.cpp``.
These conventions allow LAMMPS to automatically integrate it into the
executable when compiling and associate your new fix class with the designated
keyword when it parses the input script.
Let's write a simple fix which will print the average velocity at the end
of each timestep. First of all, implement a constructor:
.. code-block:: c++
FixPrintVel::FixPrintVel(LAMMPS *lmp, int narg, char **arg)
: Fix(lmp, narg, arg)
{
if (narg < 4)
error->all(FLERR,"Illegal fix print/vel command");
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery <= 0)
error->all(FLERR,"Illegal fix print/vel command");
}
In the constructor you should parse your fix arguments which are
specified in the script. All fixes have pretty much the same syntax:
``fix <fix-ID> <fix group> <fix name> <fix arguments ...>``. The
first 3 parameters are parsed by Fix base class constructor, while
``<fix arguments>`` should be parsed by you. In our case, we need to
specify how often we want to print an average velocity. For instance,
once in 50 timesteps: ``fix 1 print/vel 50``. There is a special variable
in the Fix class called ``nevery`` which specifies how often the method
``end_of_step()`` is called. Thus all we need to do is just set it up.
The next method we need to implement is ``setmask()``:
.. code-block:: c++
int FixPrintVel::setmask()
{
int mask = 0;
mask |= FixConst::END_OF_STEP;
return mask;
}
Here the user specifies which methods of your fix should be called
during execution. The constant ``END_OF_STEP`` corresponds to the
``end_of_step()`` method. The most important available methods that
are called during a timestep and the order in which they are called
are shown in the previous section.
.. code-block:: c++
void FixPrintVel::end_of_step()
{
// for add3, scale3
using namespace MathExtra;
double** v = atom->v;
int nlocal = atom->nlocal;
double localAvgVel[4]; // 4th element for particles count
memset(localAvgVel, 0, 4 * sizeof(double));
for (int particleInd = 0; particleInd < nlocal; ++particleInd) {
add3(localAvgVel, v[particleInd], localAvgVel);
}
localAvgVel[3] = nlocal;
double globalAvgVel[4];
memset(globalAvgVel, 0, 4 * sizeof(double));
MPI_Allreduce(localAvgVel, globalAvgVel, 4, MPI_DOUBLE, MPI_SUM, world);
scale3(1.0 / globalAvgVel[3], globalAvgVel);
if ((comm->me == 0) && screen) {
fmt::print(screen,"{}, {}, {}\n",
globalAvgVel[0], globalAvgVel[1], globalAvgVel[2]);
}
}
In the code above, we use MathExtra routines defined in
``math_extra.h``. There are bunch of math functions to work with
arrays of doubles as with math vectors. It is also important to note
that LAMMPS code should always assume to be run in parallel and that
atom data is thus distributed across the MPI ranks. Thus you can
only process data from local atoms directly and need to use MPI library
calls to combine or exchange data. For serial execution, LAMMPS
comes bundled with the MPI STUBS library that contains the MPI library
function calls in dummy versions that only work for a single MPI rank.
In this code we use an instance of Atom class. This object is stored
in the Pointers class (see ``pointers.h``) which is the base class of
the Fix base class. This object contains references to various class
instances (the original instances are created and held by the LAMMPS
class) with all global information about the simulation system.
Data from the Pointers class is available to all classes inherited from
it using protected inheritance. Hence when you write you own class,
which is going to use LAMMPS data, don't forget to inherit from Pointers
or pass an Pointer to it to all functions that need access. When writing
fixes we inherit from class Fix which is inherited from Pointers so
there is no need to inherit from it directly.
The code above computes average velocity for all particles in the
simulation. Yet you have one unused parameter in fix call from the
script: ``group_name``. This parameter specifies the group of atoms
used in the fix. So we should compute average for all particles in the
simulation only if ``group_name == "all"``, but it can be any group.
The group membership information of an atom is contained in the *mask*
property of and atom and the bit corresponding to a given group is
stored in the groupbit variable which is defined in Fix base class:
.. code-block:: c++
for (int i = 0; i < nlocal; ++i) {
if (atom->mask[i] & groupbit) {
// Do all job here
}
}
Class Atom encapsulates atoms positions, velocities, forces, etc. User
can access them using particle index. Note, that particle indexes are
usually changed every few timesteps because of neighbor list rebuilds
and spatial sorting (to improve cache efficiency).
Let us consider another Fix example: We want to have a fix which stores
atoms position from previous time step in your fix. The local atoms
indexes may not be valid on the next iteration. In order to handle
this situation there are several methods which should be implemented:
- ``double memory_usage()``: return how much memory the fix uses (optional)
- ``void grow_arrays(int)``: do reallocation of the per particle arrays in your fix
- ``void copy_arrays(int i, int j, int delflag)``: copy i-th per-particle
information to j-th. Used when atom sorting is performed. if delflag is set
and atom j owns a body, move the body information to atom i.
- ``void set_arrays(int i)``: sets i-th particle related information to zero
Note, that if your class implements these methods, it must call add calls of
add_callback and delete_callback to constructor and destructor. Since we want
to store positions of atoms from previous timestep, we need to add
``double** xold`` to the header file. Than add allocation code
to the constructor:
.. code-block:: c++
FixSavePos::FixSavePos(LAMMPS *lmp, int narg, char **arg), xold(nullptr)
{
//...
memory->create(xold, atom->nmax, 3, "FixSavePos:x");
atom->add_callback(0);
}
FixSavePos::~FixSavePos() {
atom->delete_callback(id, 0);
memory->destroy(xold);
}
Implement the aforementioned methods:
.. code-block:: c++
double FixSavePos::memory_usage()
{
int nmax = atom->nmax;
double bytes = 0.0;
bytes += nmax * 3 * sizeof(double);
return bytes;
}
void FixSavePos::grow_arrays(int nmax)
{
memory->grow(xold, nmax, 3, "FixSavePos:xold");
}
void FixSavePos::copy_arrays(int i, int j, int delflag)
{
memcpy(xold[j], xold[i], sizeof(double) * 3);
}
void FixSavePos::set_arrays(int i)
{
memset(xold[i], 0, sizeof(double) * 3);
}
int FixSavePos::pack_exchange(int i, double *buf)
{
int m = 0;
buf[m++] = xold[i][0];
buf[m++] = xold[i][1];
buf[m++] = xold[i][2];
return m;
}
int FixSavePos::unpack_exchange(int nlocal, double *buf)
{
int m = 0;
xold[nlocal][0] = buf[m++];
xold[nlocal][1] = buf[m++];
xold[nlocal][2] = buf[m++];
return m;
}
Now, a little bit about memory allocation. We use the Memory class which
is just a bunch of template functions for allocating 1D and 2D
arrays. So you need to add include ``memory.h`` to have access to them.
Finally, if you need to write/read some global information used in
your fix to the restart file, you might do it by setting flag
``restart_global = 1`` in the constructor and implementing methods void
``write_restart(FILE *fp)`` and ``void restart(char *buf)``.
If, in addition, you want to write the per-atom property to restart
files additional settings and functions are needed:
- a fix flag indicating this needs to be set ``restart_peratom = 1;``
- ``atom->add_callback()`` and ``atom->delete_callback()`` must be called
a second time with the final argument set to 1 instead of 0 (indicating
restart processing instead of per-atom data memory management).
- the functions ``void pack_restart(int i, double *buf)`` and
``void unpack_restart(int nlocal, int nth)`` need to be implemented

308
doc/src/Developer_write_pair.rst Normal file
View File

@ -0,0 +1,308 @@
Writing new pair styles
^^^^^^^^^^^^^^^^^^^^^^^
Pair styles are at the core of most simulations with LAMMPS since they
are used to compute the forces (plus energy and virial contributions, if
needed) on atoms for pairs of atoms within a given cutoff. This is
often the dominant computation in LAMMPS and sometimes even the only
one. Pair styles can be grouped in multiple categories:
#. simple pairwise additive interactions of point particles
(e.g. :doc:`Lennard-Jones <pair_lj>`, :doc:`Morse <pair_morse>`,
:doc:`Buckingham <pair_buck>`)
#. pairwise additive interactions of point particles with added
:doc:`Coulomb <pair_coul>` interactions
#. manybody interactions of point particles (e.g. :doc:`EAM <pair_eam>`,
:doc:`Tersoff <pair_tersoff>`)
#. complex interactions that include additional per-atom properties
(e.g. Discrete Element Models (DEM), Peridynamics, Ellipsoids,
Smoothed Particle Hydrodynamics (SPH))
In the text below we will discuss aspects of implementing pair styles in
LAMMPS by looking at representative case studies. The design of LAMMPS
allows to focus on the essentials, which is to compute the forces (and
energies or virial contributions), enter and manage the global settings
as well as the potential parameters, and the pair style specific parts
of reading and writing restart and data files. Most of the complex
tasks like management of the atom positions, domain decomposition and
boundaries, or neighbor list creation are handled transparently by other
parts of the LAMMPS code.
As shown on the page for :doc:`writing or extending pair styles
<Modify_pair>` for the implementation of a pair style a new class must
be written that is either directly or indirectly derived from the
``Pair`` class. In that derived class there are three *required*
methods in addition to the constructor that must be implemented since
they are "pure" in the base class: ``Pair::compute()``,
``Pair::settings()``, ``Pair::coeff()``. All other methods are optional
and have default implementations in the base class (most of which do
nothing), but they may need to be overridden depending on the
requirements of the model.
Case 1: a pairwise additive model
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In this section we will describe the complete procedure of adding
a simple pair style to LAMMPS: an empirical model that can be used
to model liquid mercury.
Model and general considerations
""""""""""""""""""""""""""""""""
The functional form of the model according to :ref:`(Bomont) <Bomont2>`
consists of a repulsive Born-Mayer exponential term and a temperature
dependent, attractive Gaussian term.
.. math::
E = A_0 \exp \left( -\alpha r \right) - A_1 \exp\left[ -\beta \left(r - r_0 \right)^2 \right]
For the application to mercury the following parameters are listed:
- :math:`A_0 = 8.2464 \times 10^{13} \; \textrm{eV}`
- :math:`\alpha = 12.48 \; \AA^{-1}`
- :math:`\beta = 0.44 \; \AA^{-2}`
- :math:`r_0 = 3.56 \; \AA`
- :math:`A_1` is temperature dependent and can be determined from
:math:`A_1 = a_0 + a_1 T + a_2 T^2` with:
- :math:`a_0 = 1.97475 \times 10^{-2} \; \textrm{eV}`
- :math:`a_1 = 8.40841 \times 10^{-5} \; \textrm{eV/K}`
- :math:`a_2 = -2.58717 \times 10^{-8} \; \textrm{eV/K}^{-2}`
With the optional cutoff, this means we have a total of 5 or 6
parameters for each pair of atom types. In addition we need to get the
temperature and a default cutoff value as global settings.
Because of the combination of Born-Mayer with a Gaussian, the pair style
shall be named "born/gauss" and thus the class name would be
``PairBornGauss`` and the source files ``pair_born_gauss.h`` and
``pair_born_gauss.cpp``. Since this is a rather uncommon potential, it
shall be added to the :ref:`EXTRA-PAIR <PKG-EXTRA-PAIR>` package. For
the implementation, we will use :doc:`pair style morse <pair_morse>` as
template.
Header file
"""""""""""
The first segment of any LAMMPS source should be the copyright and
license statement. Note the marker in the first line to indicate to
editors like emacs that this file is a C++ source, even though the .h
extension suggests a C source.
.. code-block:: c++
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
Every pair style must be registered in LAMMPS by writing the following
lines of code in the second part of the header before the include guards
for the class definition:
.. code-block:: c++
#ifdef PAIR_CLASS
// clang-format off
PairStyle(born/gauss,PairBornGauss);
// clang-format on
#else
/* the definition of the PairBornGauss class (see below) is inserted here */
#endif
This second segment of the header file will be included by the ``Force``
class in ``force.cpp`` to build a map of "factory functions" that will
create an instance of these classes and return a pointer to it. The map
connects the name of the pair style, "born/gauss", to the name of the
class, ``PairBornGauss``. The list of header files to include is
automatically updated by the build system, so the presence of the file
in the ``src/EXTRA-PAIR`` folder and the enabling of the EXTRA-PAIR
package will trigger that LAMMPS includes the new pair style when it is
(re-)compiled. The "// clang-format" format comments are needed so that
running :ref:`clang-format <clang-format>` on the file will not insert
blanks between "born", "/", and "gauss" which would break the
``PairStyle`` macro.
The third segment of the header is the actual class definition of the
``PairBornGauss`` class. This has the prototypes for all member
functions that will be implemented by this pair style. This includes a
number of optional functions. All functions that were labeled in the
base class as "virtual" must be given the "override" property as it is
done in the code shown below. This helps to detect unexpected
mismatches as compile errors in case the signature of a function is
changed in the base class. For example, if this change would add an
optional argument with a default value, then all existing source code
calling the function would not need changes and still compile, but the
function in the derived class would no longer override the one in the
base class due to the different number of arguments and the behavior of
the pair style is thus changed in an unintended way. Using "override"
prevents such issues.
Also variables and arrays for storing global settings and potential
parameters are defined. Since those are internal to the class, they are
placed after a "protected:" label.
.. code-block:: c++
#ifndef LMP_PAIR_BORN_GAUSS_H
#define LMP_PAIR_BORN_GAUSS_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBornGauss : public Pair {
public:
PairBornGauss(class LAMMPS *);
~PairBornGauss() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
void write_data(FILE *) override;
void write_data_all(FILE *) override;
double single(int, int, int, int, double, double, double, double &) override;
void *extract(const char *, int &) override;
protected:
double cut_global, temperature;
double **cut;
double **biga0, **alpha, **biga1, **beta, **r0;
double **a0, **a1, **a2;
double **offset;
virtual void allocate();
};
} // namespace LAMMPS_NS
#endif
Some details of the class definition will be discussed later.
Implementation file
"""""""""""""""""""
We move on to the implementation in the ``pair_born_gauss.cpp`` file.
This file also starts with the LAMMPS copyright and license header.
Below is typically the space where comments may be added with additional
information about this specific file, the author(s) and affiliation(s)
and email address(es) so the author can be easily contacted in case
there are questions about the implementation later. Since the file(s)
may be around for a long time, it is beneficial to use some kind of
"permanent" email address, if possible.
.. code-block:: c++
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
// Contributing author: Axel Kohlmeyer, Temple University, akohlmey@gmail.com
#include "pair_born_gauss.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
The second section of the implementation file has various include
statements. The include file for the class header has to come first,
then LAMMPS classes (sorted alphabetically) and system headers and
others, if needed. Note the standardized C++ notation for headers of
C-library functions. The final statement of this segment imports the
``LAMMPS_NS::`` namespace globally for this file. This way, all LAMMPS
specific functions and classes do not have to be prefixed with
``LAMMPS_NS::``.
Constructor and destructor (required)
"""""""""""""""""""""""""""""""""""""
The first two functions in the implementation source file are typically
the constructor and the destructor.
.. code-block:: c++
/* ---------------------------------------------------------------------- */
PairBornGauss::PairBornGauss(LAMMPS *lmp) : Pair(lmp), temp(nullptr)
{
single_enable = 1;
respa_enable = 0;
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairBornGauss::~PairBornGauss()
{
}
In the constructor you should parse your fix arguments which are
specified in the script. All fixes have pretty much the same syntax:
``fix <fix-ID> <fix group> <fix name> <fix arguments ...>``. The
first 3 parameters are parsed by Fix base class constructor, while
``<fix arguments>`` should be parsed by you. In our case, we need to
specify how often we want to print an average velocity. For instance,
once in 50 timesteps: ``fix 1 print/vel 50``. There is a special variable
in the Fix class called ``nevery`` which specifies how often the method
``end_of_step()`` is called. Thus all we need to do is just set it up.
Settings and coefficients (required)
""""""""""""""""""""""""""""""""""""
The next method we need to implement is ``setmask()``:
.. code-block:: c++
int FixPrintVel::setmask()
{
int mask = 0;
mask |= FixConst::END_OF_STEP;
return mask;
}
Computing forces from the neighbor list (required)
""""""""""""""""""""""""""""""""""""""""""""""""""
xxxxxadfasdf
--------------
.. _Bomont2:
**(Bomont)** Bomont, Bretonnet, J. Chem. Phys. 124, 054504 (2006)