Files
lammps/doc/html/Developer_write_pair.html
2025-01-13 14:55:48 +00:00

1438 lines
203 KiB
HTML
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

<!DOCTYPE html>
<html class="writer-html5" lang="en" >
<head>
<meta charset="utf-8" /><meta name="viewport" content="width=device-width, initial-scale=1" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<title>4.8.1. Writing new pair styles &mdash; LAMMPS documentation</title>
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
<link rel="stylesheet" href="_static/css/theme.css" type="text/css" />
<link rel="stylesheet" href="_static/sphinx-design.min.css" type="text/css" />
<link rel="stylesheet" href="_static/css/lammps.css" type="text/css" />
<link rel="shortcut icon" href="_static/lammps.ico"/>
<link rel="canonical" href="https://docs.lammps.org/Developer_write_pair.html" />
<!--[if lt IE 9]>
<script src="_static/js/html5shiv.min.js"></script>
<![endif]-->
<script src="_static/jquery.js?v=5d32c60e"></script>
<script src="_static/_sphinx_javascript_frameworks_compat.js?v=2cd50e6c"></script>
<script src="_static/documentation_options.js?v=5929fcd5"></script>
<script src="_static/doctools.js?v=9bcbadda"></script>
<script src="_static/sphinx_highlight.js?v=dc90522c"></script>
<script src="_static/design-tabs.js?v=f930bc37"></script>
<script async="async" src="_static/mathjax/es5/tex-mml-chtml.js?v=cadf963e"></script>
<script src="_static/js/theme.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
<link rel="next" title="4.8.7. Writing a new fix style" href="Developer_write_fix.html" />
<link rel="prev" title="4.8. Writing new styles" href="Developer_write.html" />
</head>
<body class="wy-body-for-nav">
<div class="wy-grid-for-nav">
<nav data-toggle="wy-nav-shift" class="wy-nav-side">
<div class="wy-side-scroll">
<div class="wy-side-nav-search" >
<a href="Manual.html">
<img src="_static/lammps-logo.png" class="logo" alt="Logo"/>
</a>
<div class="lammps_version">Version: <b>19 Nov 2024</b></div>
<div class="lammps_release">git info: </div>
<div role="search">
<form id="rtd-search-form" class="wy-form" action="search.html" method="get">
<input type="text" name="q" placeholder="Search docs" aria-label="Search docs" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
</div>
</div><div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="Navigation menu">
<p class="caption" role="heading"><span class="caption-text">User Guide</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Intro.html">1. Introduction</a></li>
<li class="toctree-l1"><a class="reference internal" href="Install.html">2. Install LAMMPS</a></li>
<li class="toctree-l1"><a class="reference internal" href="Build.html">3. Build LAMMPS</a></li>
<li class="toctree-l1"><a class="reference internal" href="Run_head.html">4. Run LAMMPS</a></li>
<li class="toctree-l1"><a class="reference internal" href="Commands.html">5. Commands</a></li>
<li class="toctree-l1"><a class="reference internal" href="Packages.html">6. Optional packages</a></li>
<li class="toctree-l1"><a class="reference internal" href="Speed.html">7. Accelerate performance</a></li>
<li class="toctree-l1"><a class="reference internal" href="Howto.html">8. Howto discussions</a></li>
<li class="toctree-l1"><a class="reference internal" href="Examples.html">9. Example scripts</a></li>
<li class="toctree-l1"><a class="reference internal" href="Tools.html">10. Auxiliary tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="Errors.html">11. Errors</a></li>
</ul>
<p class="caption" role="heading"><span class="caption-text">Programmer Guide</span></p>
<ul class="current">
<li class="toctree-l1"><a class="reference internal" href="Library.html">1. LAMMPS Library Interfaces</a></li>
<li class="toctree-l1"><a class="reference internal" href="Python_head.html">2. Use Python with LAMMPS</a></li>
<li class="toctree-l1"><a class="reference internal" href="Modify.html">3. Modifying &amp; extending LAMMPS</a></li>
<li class="toctree-l1 current"><a class="reference internal" href="Developer.html">4. Information for Developers</a><ul class="current">
<li class="toctree-l2"><a class="reference internal" href="Developer_org.html">4.1. Source files</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_org.html#class-topology">4.2. Class topology</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_code_design.html">4.3. Code design</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_parallel.html">4.4. Parallel algorithms</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_atom.html">4.5. Accessing per-atom data</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_comm_ops.html">4.6. Communication patterns</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_flow.html">4.7. How a timestep works</a></li>
<li class="toctree-l2 current"><a class="reference internal" href="Developer_write.html">4.8. Writing new styles</a><ul class="current">
<li class="toctree-l3 current"><a class="current reference internal" href="#">4.8.1. Writing new pair styles</a></li>
<li class="toctree-l3"><a class="reference internal" href="#package-and-build-system-considerations">4.8.2. Package and build system considerations</a></li>
<li class="toctree-l3"><a class="reference internal" href="#case-1-a-pairwise-additive-model">4.8.3. Case 1: a pairwise additive model</a></li>
<li class="toctree-l3"><a class="reference internal" href="#case-2-a-many-body-potential">4.8.4. Case 2: a many-body potential</a></li>
<li class="toctree-l3"><a class="reference internal" href="#case-3-a-potential-requiring-communication">4.8.5. Case 3: a potential requiring communication</a></li>
<li class="toctree-l3"><a class="reference internal" href="#case-4-potentials-without-a-compute-function">4.8.6. Case 4: potentials without a compute() function</a></li>
<li class="toctree-l3"><a class="reference internal" href="Developer_write_fix.html">4.8.7. Writing a new fix style</a></li>
<li class="toctree-l3"><a class="reference internal" href="Developer_write_command.html">4.8.8. Writing a new command style</a></li>
<li class="toctree-l3"><a class="reference internal" href="Developer_write_command.html#case-1-implementing-the-geturl-command">4.8.9. Case 1: Implementing the geturl command</a></li>
</ul>
</li>
<li class="toctree-l2"><a class="reference internal" href="Developer_notes.html">4.9. Notes for developers and code maintainers</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_updating.html">4.10. Notes for updating code written for older LAMMPS versions</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_plugins.html">4.11. Writing plugins</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_unittest.html">4.12. Adding tests for unit testing</a></li>
<li class="toctree-l2"><a class="reference internal" href="Classes.html">4.13. C++ base classes</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_platform.html">4.14. Platform abstraction functions</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html">4.15. Utility functions</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html#special-math-functions">4.16. Special Math functions</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html#tokenizer-classes">4.17. Tokenizer classes</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html#argument-parsing-classes">4.18. Argument parsing classes</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html#file-reader-classes">4.19. File reader classes</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html#memory-pool-classes">4.20. Memory pool classes</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html#eigensolver-functions">4.21. Eigensolver functions</a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_utils.html#communication-buffer-coding-with-ubuf">4.22. Communication buffer coding with <em>ubuf</em></a></li>
<li class="toctree-l2"><a class="reference internal" href="Developer_grid.html">4.23. Use of distributed grids within style classes</a></li>
</ul>
</li>
</ul>
<p class="caption" role="heading"><span class="caption-text">Command Reference</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="commands_list.html">Commands</a></li>
<li class="toctree-l1"><a class="reference internal" href="fixes.html">Fix Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="computes.html">Compute Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="pairs.html">Pair Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="bonds.html">Bond Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="angles.html">Angle Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="dihedrals.html">Dihedral Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="impropers.html">Improper Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="dumps.html">Dump Styles</a></li>
<li class="toctree-l1"><a class="reference internal" href="fix_modify_atc_commands.html">fix_modify AtC commands</a></li>
<li class="toctree-l1"><a class="reference internal" href="Bibliography.html">Bibliography</a></li>
</ul>
</div>
</div>
</nav>
<section data-toggle="wy-nav-shift" class="wy-nav-content-wrap"><nav class="wy-nav-top" aria-label="Mobile navigation menu" >
<i data-toggle="wy-nav-top" class="fa fa-bars"></i>
<a href="Manual.html">LAMMPS</a>
</nav>
<div class="wy-nav-content">
<div class="rst-content style-external-links">
<div role="navigation" aria-label="Page navigation">
<ul class="wy-breadcrumbs">
<li><a href="Manual.html" class="icon icon-home" aria-label="Home"></a></li>
<li class="breadcrumb-item"><a href="Developer.html"><span class="section-number">4. </span>Information for Developers</a></li>
<li class="breadcrumb-item"><a href="Developer_write.html"><span class="section-number">4.8. </span>Writing new styles</a></li>
<li class="breadcrumb-item active"><span class="section-number">4.8.1. </span>Writing new pair styles</li>
<li class="wy-breadcrumbs-aside">
<a href="https://www.lammps.org"><img src="_static/lammps-logo.png" width="64" height="16" alt="LAMMPS Homepage"></a> | <a href="Commands_all.html">Commands</a>
</li>
</ul><div class="rst-breadcrumbs-buttons" role="navigation" aria-label="Sequential page navigation">
<a href="Developer_write.html" class="btn btn-neutral float-left" title="4.8. Writing new styles" accesskey="p"><span class="fa fa-arrow-circle-left" aria-hidden="true"></span> Previous</a>
<a href="Developer_write_fix.html" class="btn btn-neutral float-right" title="4.8.7. Writing a new fix style" accesskey="n">Next <span class="fa fa-arrow-circle-right" aria-hidden="true"></span></a>
</div>
<hr/>
</div>
<div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
<div itemprop="articleBody">
<p><span class="math notranslate nohighlight">\(\renewcommand{\AA}{\text{Å}}\)</span></p>
<section id="writing-new-pair-styles">
<h1><span class="section-number">4.8.1. </span>Writing new pair styles<a class="headerlink" href="#writing-new-pair-styles" title="Link to this heading"></a></h1>
<p>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 or small clusters 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 into multiple categories:</p>
<ol class="arabic simple">
<li><p>simple pairwise additive interactions of point particles
(e.g. <a class="reference internal" href="pair_lj.html"><span class="doc">Lennard-Jones</span></a>, <a class="reference internal" href="pair_morse.html"><span class="doc">Morse</span></a>,
<a class="reference internal" href="pair_buck.html"><span class="doc">Buckingham</span></a>)</p></li>
<li><p>pairwise additive interactions of point particles with added
<a class="reference internal" href="pair_coul.html"><span class="doc">Coulomb</span></a> interactions or <em>only</em> the Coulomb interactions</p></li>
<li><p>manybody interactions of point particles (e.g. <a class="reference internal" href="pair_eam.html"><span class="doc">EAM</span></a>,
<a class="reference internal" href="pair_tersoff.html"><span class="doc">Tersoff</span></a>)</p></li>
<li><p>complex interactions that include additional per-atom properties
(e.g. Discrete Element Models (DEM), Peridynamics, Ellipsoids)</p></li>
<li><p>special purpose pair styles that may not even compute forces like
<a class="reference internal" href="pair_zero.html"><span class="doc">pair_style zero</span></a> and <a class="reference internal" href="pair_tracker.html"><span class="doc">pair_style tracker</span></a>, or are a wrapper for multiple kinds of interactions
like <a class="reference internal" href="pair_hybrid.html"><span class="doc">pair_style hybrid</span></a>, <a class="reference internal" href="pair_list.html"><span class="doc">pair_style list</span></a>,
and <a class="reference internal" href="pair_kim.html"><span class="doc">pair_style kim</span></a></p></li>
</ol>
<p>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 developers 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.</p>
<p>As shown on the page for <a class="reference internal" href="Modify_pair.html"><span class="doc">writing or extending pair styles</span></a>, in order to implement a new pair style, a new class must
be written that is either directly or indirectly derived from the
<code class="docutils literal notranslate"><span class="pre">Pair</span></code> class. If that class is directly derived from <code class="docutils literal notranslate"><span class="pre">Pair</span></code>, there
are three methods that <em>must</em> be re-implemented, since they are “pure”
in the base class: <code class="docutils literal notranslate"><span class="pre">Pair::compute()</span></code>, <code class="docutils literal notranslate"><span class="pre">Pair::settings()</span></code>,
<code class="docutils literal notranslate"><span class="pre">Pair::coeff()</span></code>. In addition, a custom constructor is needed. 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.</p>
<p>We are looking at the following cases:</p>
<ul class="simple">
<li><p><a class="reference internal" href="#case-1-a-pairwise-additive-model">Case 1: a pairwise additive model</a></p></li>
<li><p><a class="reference internal" href="#case-2-a-many-body-potential">Case 2: a many-body potential</a></p></li>
<li><p><a class="reference internal" href="#case-3-a-potential-requiring-communication">Case 3: a potential requiring communication</a></p></li>
<li><p><a class="reference internal" href="#case-4-potentials-without-a-compute-function">Case 4: potentials without a compute() function</a></p></li>
</ul>
</section>
<section id="package-and-build-system-considerations">
<h1><span class="section-number">4.8.2. </span>Package and build system considerations<a class="headerlink" href="#package-and-build-system-considerations" title="Link to this heading"></a></h1>
<p>In general, new pair styles should be added to the <a class="reference internal" href="Packages_details.html#pkg-extra-pair"><span class="std std-ref">EXTRA-PAIR
package</span></a> unless they are an accelerated pair style and
then they should be added to the corresponding accelerator package
(<a class="reference internal" href="Packages_details.html#pkg-gpu"><span class="std std-ref">GPU</span></a>, <a class="reference internal" href="Packages_details.html#pkg-intel"><span class="std std-ref">INTEL</span></a>, <a class="reference internal" href="Packages_details.html#pkg-kokkos"><span class="std std-ref">KOKKOS</span></a>, <a class="reference internal" href="Packages_details.html#pkg-openmp"><span class="std std-ref">OPENMP</span></a>, <a class="reference internal" href="Packages_details.html#pkg-opt"><span class="std std-ref">OPT</span></a>). If
you feel that your contribution should be added to a different package,
please consult with the LAMMPS developers first.</p>
<p>The contributed code needs to support the <a class="reference internal" href="Build_make.html"><span class="doc">traditional GNU make
build process</span></a> <strong>and</strong> the <a class="reference internal" href="Build_cmake.html"><span class="doc">CMake build process</span></a>. For the GNU make process and if the package has an
<code class="docutils literal notranslate"><span class="pre">Install.sh</span></code> file, most likely that file needs to be updated to
correctly copy the sources when installing the package and properly
delete them when uninstalling. This is particularly important when
added a new pair style that is a derived class from an existing pair
style in a package, so that its installation depends on the the
installation status of the package of the derived class. For the CMake
process, it is sometimes necessary to make changes to the package
specific CMake scripting in <code class="docutils literal notranslate"><span class="pre">cmake/Modules/Packages</span></code>.</p>
</section>
<hr class="docutils" />
<section id="case-1-a-pairwise-additive-model">
<h1><span class="section-number">4.8.3. </span>Case 1: a pairwise additive model<a class="headerlink" href="#case-1-a-pairwise-additive-model" title="Link to this heading"></a></h1>
<p>In this section, we will describe the procedure of adding a simple pair
style to LAMMPS: an empirical model that can be used to model liquid
mercury. The pair style shall be called <a class="reference internal" href="pair_born_gauss.html"><span class="doc">bond/gauss</span></a> and the complete implementation can be found in the
files <code class="docutils literal notranslate"><span class="pre">src/EXTRA-PAIR/pair_born_gauss.cpp</span></code> and
<code class="docutils literal notranslate"><span class="pre">src/EXTRA-PAIR/pair_born_gauss.h</span></code> of the LAMMPS source code.</p>
<section id="model-and-general-considerations">
<h2>Model and general considerations<a class="headerlink" href="#model-and-general-considerations" title="Link to this heading"></a></h2>
<p>The functional form of the model according to <a class="reference internal" href="#bomont2"><span class="std std-ref">(Bomont)</span></a>
consists of a repulsive Born-Mayer exponential term and a temperature
dependent, attractive Gaussian term.</p>
<div class="math notranslate nohighlight">
\[E = A_0 \exp \left( -\alpha r \right) - A_1 \exp\left[ -\beta \left(r - r_0 \right)^2 \right]\]</div>
<p>For the application to mercury, the following parameters are listed:</p>
<ul class="simple">
<li><p><span class="math notranslate nohighlight">\(A_0 = 8.2464 \times 10^{13} \; \textrm{eV}\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(\alpha = 12.48 \; \AA^{-1}\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(\beta = 0.44 \; \AA^{-2}\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(r_0 = 3.56 \; \AA\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(A_1\)</span> is temperature dependent and can be determined from
<span class="math notranslate nohighlight">\(A_1 = a_0 + a_1 T + a_2 T^2\)</span> with:</p>
<ul>
<li><p><span class="math notranslate nohighlight">\(a_0 = 1.97475 \times 10^{-2} \; \textrm{eV}\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(a_1 = 8.40841 \times 10^{-5} \; \textrm{eV/K}\)</span></p></li>
<li><p><span class="math notranslate nohighlight">\(a_2 = -2.58717 \times 10^{-8} \; \textrm{eV/K}^{-2}\)</span></p></li>
</ul>
</li>
</ul>
<p>With the optional cutoff, this means we have a total of 5 or 6
parameters for each pair of atom types. Additionally, we need to input a
default cutoff value as a global setting.</p>
<p>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
<code class="docutils literal notranslate"><span class="pre">PairBornGauss</span></code> and the source files <code class="docutils literal notranslate"><span class="pre">pair_born_gauss.h</span></code> and
<code class="docutils literal notranslate"><span class="pre">pair_born_gauss.cpp</span></code>. Since this is a rather uncommon potential, it
shall be added to the <a class="reference internal" href="Packages_details.html#pkg-extra-pair"><span class="std std-ref">EXTRA-PAIR</span></a> package.</p>
</section>
<section id="header-file">
<h2>Header file<a class="headerlink" href="#header-file" title="Link to this heading"></a></h2>
<p>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 (this is a convention inherited from the
very beginning of the C++ version of LAMMPS).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* -*- c++ -*- ----------------------------------------------------------</span>
<span class="cm"> LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator</span>
<span class="cm"> https://www.lammps.org/, Sandia National Laboratories</span>
<span class="cm"> LAMMPS development team: developers@lammps.org</span>
<span class="cm"> Copyright (2003) Sandia Corporation. Under the terms of Contract</span>
<span class="cm"> DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains</span>
<span class="cm"> certain rights in this software. This software is distributed under</span>
<span class="cm"> the GNU General Public License.</span>
<span class="cm"> See the README file in the top-level LAMMPS directory.</span>
<span class="cm">------------------------------------------------------------------------- */</span>
</pre></div>
</div>
<p>Every pair style must be registered in LAMMPS by including the following
lines of code in the second part of the header after the copyright
message and before the include guards for the class definition:</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cp">#ifdef PAIR_CLASS</span>
<span class="c1">// clang-format off</span>
<span class="n">PairStyle</span><span class="p">(</span><span class="n">born</span><span class="o">/</span><span class="n">gauss</span><span class="p">,</span><span class="n">PairBornGauss</span><span class="p">);</span>
<span class="c1">// clang-format on</span>
<span class="cp">#else</span>
<span class="cm">/* the definition of the PairBornGauss class (see below) is inserted here */</span>
<span class="cp">#endif</span>
</pre></div>
</div>
<p>This block between <code class="docutils literal notranslate"><span class="pre">#ifdef</span> <span class="pre">PAIR_CLASS</span></code> and <code class="docutils literal notranslate"><span class="pre">#else</span></code> will be
included by the <code class="docutils literal notranslate"><span class="pre">Force</span></code> class in <code class="docutils literal notranslate"><span class="pre">force.cpp</span></code> 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, <code class="docutils literal notranslate"><span class="pre">PairBornGauss</span></code>. During
compilation, LAMMPS constructs a file <code class="docutils literal notranslate"><span class="pre">style_pair.h</span></code> that contains
<code class="docutils literal notranslate"><span class="pre">#include</span></code> statements for all “installed” pair styles. Before
including <code class="docutils literal notranslate"><span class="pre">style_pair.h</span></code> into <code class="docutils literal notranslate"><span class="pre">force.cpp</span></code>, the <code class="docutils literal notranslate"><span class="pre">PAIR_CLASS</span></code> define
is set and the <code class="docutils literal notranslate"><span class="pre">PairStyle(name,class)</span></code> macro defined. The code of the
macro adds the installed pair styles to the “factory map” which enables
the <a class="reference internal" href="pair_style.html"><span class="doc">pair_style command</span></a> to create the pair style
instance.</p>
<p>The list of header files to include is automatically updated by the
build system if there are new files, so the presence of the new header
file in the <code class="docutils literal notranslate"><span class="pre">src/EXTRA-PAIR</span></code> folder and the enabling of the EXTRA-PAIR
package will trigger LAMMPS to include the new pair style when it is
(re-)compiled. The “// clang-format” format comments are needed so that
running <a class="reference internal" href="Build_development.html#clang-format"><span class="std std-ref">clang-format</span></a> on the file will not insert
unwanted blanks between “born”, “/”, and “gauss” which would break the
<code class="docutils literal notranslate"><span class="pre">PairStyle</span></code> macro.</p>
<p>The third part of the header file is the actual class definition of the
<code class="docutils literal notranslate"><span class="pre">PairBornGauss</span></code> class. This has the prototypes for all member
functions that will be implemented by this pair style. This includes
<a class="reference internal" href="Modify_pair.html"><span class="doc">a few required and a number of optional functions</span></a>.
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.</p>
<p>The “override” property helps to detect unexpected mismatches because
compilation will stop with an error in case the signature of a function
is changed in the base class without also changing it in all derived
classes. For example, if this change added an optional argument with a
default value, then all existing source code <em>calling</em> 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 the “override” keyword
prevents such issues.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cp">#ifndef LMP_PAIR_BORN_GAUSS_H</span>
<span class="cp">#define LMP_PAIR_BORN_GAUSS_H</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;pair.h&quot;</span>
<span class="k">namespace</span><span class="w"> </span><span class="nn">LAMMPS_NS</span><span class="w"> </span><span class="p">{</span>
<span class="k">class</span><span class="w"> </span><span class="nc">PairBornGauss</span><span class="w"> </span><span class="o">:</span><span class="w"> </span><span class="k">public</span><span class="w"> </span><span class="n">Pair</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">public</span><span class="o">:</span>
<span class="w"> </span><span class="n">PairBornGauss</span><span class="p">(</span><span class="k">class</span><span class="w"> </span><span class="nc">LAMMPS</span><span class="w"> </span><span class="o">*</span><span class="p">);</span>
<span class="w"> </span><span class="o">~</span><span class="n">PairBornGauss</span><span class="p">()</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">compute</span><span class="p">(</span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">settings</span><span class="p">(</span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">coeff</span><span class="p">(</span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="nf">init_one</span><span class="p">(</span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">write_restart</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">read_restart</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">write_restart_settings</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">read_restart_settings</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">write_data</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">write_data_all</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="nf">single</span><span class="p">(</span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">&amp;</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="o">*</span><span class="nf">extract</span><span class="p">(</span><span class="k">const</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">*</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="o">&amp;</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
</pre></div>
</div>
<p>Also, variables and arrays for storing global settings and potential
parameters are defined. Since these are internal to the class, they are
placed after a “protected:” label.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="w"> </span><span class="k">protected</span><span class="o">:</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">cut_global</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">**</span><span class="n">cut</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">**</span><span class="n">biga0</span><span class="p">,</span><span class="w"> </span><span class="o">**</span><span class="n">alpha</span><span class="p">,</span><span class="w"> </span><span class="o">**</span><span class="n">biga1</span><span class="p">,</span><span class="w"> </span><span class="o">**</span><span class="n">beta</span><span class="p">,</span><span class="w"> </span><span class="o">**</span><span class="n">r0</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">**</span><span class="n">a0</span><span class="p">,</span><span class="w"> </span><span class="o">**</span><span class="n">a1</span><span class="p">,</span><span class="w"> </span><span class="o">**</span><span class="n">a2</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">**</span><span class="n">offset</span><span class="p">;</span>
<span class="w"> </span><span class="k">virtual</span><span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">allocate</span><span class="p">();</span>
<span class="p">};</span>
<span class="p">}</span><span class="w"> </span><span class="c1">// namespace LAMMPS_NS</span>
<span class="cp">#endif</span>
</pre></div>
</div>
</section>
<section id="implementation-file">
<h2>Implementation file<a class="headerlink" href="#implementation-file" title="Link to this heading"></a></h2>
<p>We move on to the implementation of the <code class="docutils literal notranslate"><span class="pre">PairBornGauss</span></code> class in the
<code class="docutils literal notranslate"><span class="pre">pair_born_gauss.cpp</span></code> file. This file also starts with a LAMMPS
copyright and license header. Below that notice is typically the space
where comments may be added with additional information about this
specific file, the author(s), affiliation(s), and email address(es).
This way the contributing author(s) can be easily contacted, when
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.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator</span>
<span class="cm"> https://www.lammps.org/, Sandia National Laboratories</span>
<span class="cm"> LAMMPS development team: developers@lammps.org</span>
<span class="cm"> Copyright (2003) Sandia Corporation. Under the terms of Contract</span>
<span class="cm"> DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains</span>
<span class="cm"> certain rights in this software. This software is distributed under</span>
<span class="cm"> the GNU General Public License.</span>
<span class="cm"> See the README file in the top-level LAMMPS directory.</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="c1">// Contributing author: Axel Kohlmeyer, Temple University, akohlmey@gmail.com</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;pair_born_gauss.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;atom.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;comm.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;error.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;fix.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;force.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;memory.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&quot;neigh_list.h&quot;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&lt;cmath&gt;</span>
<span class="cp">#include</span><span class="w"> </span><span class="cpf">&lt;cstring&gt;</span>
<span class="k">using</span><span class="w"> </span><span class="k">namespace</span><span class="w"> </span><span class="nn">LAMMPS_NS</span><span class="p">;</span>
</pre></div>
</div>
<p>The second section of the implementation file has various include
statements. The include file for the class header has to come first,
then a block of LAMMPS classes (sorted alphabetically) followed by a
block of system headers and others, if needed. Note the standardized
C++ notation for headers of C-library functions (<code class="docutils literal notranslate"><span class="pre">cmath</span></code> instead of
<code class="docutils literal notranslate"><span class="pre">math.h</span></code>). The final statement of this segment imports the
<code class="docutils literal notranslate"><span class="pre">LAMMPS_NS::</span></code> namespace globally for this file. This way, all LAMMPS
specific functions and classes do not have to be prefixed with
<code class="docutils literal notranslate"><span class="pre">LAMMPS_NS::</span></code>.</p>
</section>
<section id="constructor-and-destructor-required">
<h2>Constructor and destructor (required)<a class="headerlink" href="#constructor-and-destructor-required" title="Link to this heading"></a></h2>
<p>The first two functions in the implementation source file are typically
the constructor and the destructor.</p>
<p>Pair styles are different from most classes in LAMMPS that define a
“style”, as their constructor only uses the LAMMPS class instance
pointer as an argument, but <strong>not</strong> the arguments of the
<a class="reference internal" href="pair_style.html"><span class="doc">pair_style command</span></a>. Instead, those arguments are
processed in the <code class="docutils literal notranslate"><span class="pre">Pair::settings()</span></code> function (or rather the version in
the derived class). The constructor is the place where global defaults
are set and specifically flags are set indicating which optional
features of a pair style are available.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ---------------------------------------------------------------------- */</span>
<span class="n">PairBornGauss</span><span class="o">::</span><span class="n">PairBornGauss</span><span class="p">(</span><span class="n">LAMMPS</span><span class="w"> </span><span class="o">*</span><span class="n">lmp</span><span class="p">)</span><span class="w"> </span><span class="o">:</span><span class="w"> </span><span class="n">Pair</span><span class="p">(</span><span class="n">lmp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">writedata</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="p">}</span>
</pre></div>
</div>
<p>The <cite>writedata = 1;</cite> statement indicates that the pair style is capable
of writing the current pair coefficient parameters to data files. That
is, the class implements specific versions for <code class="docutils literal notranslate"><span class="pre">Pair::data()</span></code> and
<code class="docutils literal notranslate"><span class="pre">Pair::data_all()</span></code>. Other statements that could be added here would
be <cite>single_enable = 1;</cite> or <cite>respa_enable = 0;</cite> to indicate that the
<code class="docutils literal notranslate"><span class="pre">Pair::single()</span></code> function is present and the
<code class="docutils literal notranslate"><span class="pre">Pair::compute_(inner|middle|outer)</span></code> functions are not, but those are
also the default settings and already set in the base class.</p>
<p>In the destructor, we need to delete all memory that was allocated by the
pair style, usually to hold force field parameters that were entered
with the <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff command</span></a>. Most of those array
pointers will need to be declared in the derived class header, but some
(e.g. setflag, cutsq) are already declared in the base class.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="n">PairBornGauss</span><span class="o">::~</span><span class="n">PairBornGauss</span><span class="p">()</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">allocated</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">setflag</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">cutsq</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">cut</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">biga0</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">alpha</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">biga1</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">beta</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">r0</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">offset</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="settings-and-coefficients-required">
<h2>Settings and coefficients (required)<a class="headerlink" href="#settings-and-coefficients-required" title="Link to this heading"></a></h2>
<p>To enter the global pair style settings and the pair style parameters,
the functions <code class="docutils literal notranslate"><span class="pre">Pair::settings()</span></code> and <code class="docutils literal notranslate"><span class="pre">Pair::coeff()</span></code> need to be
re-implemented. The arguments to the <code class="docutils literal notranslate"><span class="pre">settings()</span></code> function are the
arguments given to the <a class="reference internal" href="pair_style.html"><span class="doc">pair_style command</span></a>.
Normally, those would already be processed as part of the constructor,
but moving this to a separate function allows users to change global
settings like the default cutoff without having to reissue all
pair_coeff commands or re-read the <code class="docutils literal notranslate"><span class="pre">Pair</span> <span class="pre">Coeffs</span></code> sections from the
data file. In the <code class="docutils literal notranslate"><span class="pre">settings()</span></code> function, also the arrays for storing
parameters, to define cutoffs, track which pairs of parameters have been
explicitly set and allocated and, if needed, initialized. In this case,
the memory allocation and initialization are moved to a function
<code class="docutils literal notranslate"><span class="pre">allocate()</span></code>.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> allocate all arrays</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::allocate</span><span class="p">()</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">allocated</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">np1</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">setflag</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:setflag&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">np1</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">i</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">np1</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">cutsq</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:cutsq&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">cut</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:cut&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">biga0</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:biga0&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">alpha</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:alpha&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">biga1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:biga1&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">beta</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:beta&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">r0</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:r0&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">offset</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="n">np1</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;pair:offset&quot;</span><span class="p">);</span>
<span class="p">}</span>
<span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> global settings</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::settings</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">narg</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="n">arg</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">narg</span><span class="w"> </span><span class="o">!=</span><span class="w"> </span><span class="mi">1</span><span class="p">)</span><span class="w"> </span><span class="n">error</span><span class="o">-&gt;</span><span class="n">all</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;Pair style bond/gauss must have exactly one argument&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">cut_global</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">numeric</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="w"> </span><span class="nb">false</span><span class="p">,</span><span class="w"> </span><span class="n">lmp</span><span class="p">);</span>
<span class="w"> </span><span class="c1">// reset per-type pair cutoffs that have been explicitly set previously</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">allocated</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">i</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="o">++</span><span class="p">)</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">])</span><span class="w"> </span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">cut_global</span><span class="p">;</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">}</span>
</pre></div>
</div>
<p>The arguments to the <code class="docutils literal notranslate"><span class="pre">coeff()</span></code> function are the arguments to the
<a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff command</span></a>. The function is also called
when processing the <code class="docutils literal notranslate"><span class="pre">Pair</span> <span class="pre">Coeffs</span></code> or <code class="docutils literal notranslate"><span class="pre">PairIJ</span> <span class="pre">Coeffs</span></code> sections of
data files. In the case of the <code class="docutils literal notranslate"><span class="pre">Pair</span> <span class="pre">Coeffs</span></code> section, there is only
one atom type per line and thus the first argument is duplicated. Since
the atom type arguments of the <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff command</span></a>
may be a range (e.g. *3 for atom types 1, 2, and 3), the
corresponding arguments are passed to the <a class="reference internal" href="Developer_utils.html#_CPPv4I0EN9LAMMPS_NS5utils6boundsEvPKciRKNSt6stringE6bigint6bigintR4TYPER4TYPEP5Error" title="LAMMPS_NS::utils::bounds"><code class="xref cpp cpp-func docutils literal notranslate"><span class="pre">utils::bounds()</span></code></a> function which will then return the low
and high end of the range. Note that the <code class="docutils literal notranslate"><span class="pre">setflag</span></code> array is set to 1
for all pairs of atom types processed by this call. This information is
later used in the <code class="docutils literal notranslate"><span class="pre">init_one()</span></code> function to determine if any coefficients
are missing and, if supported by the potential, generate those missing
coefficients from the selected mixing rule.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> set coeffs for one or more type pairs</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::coeff</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">narg</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="n">arg</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">narg</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="mi">7</span><span class="w"> </span><span class="o">||</span><span class="w"> </span><span class="n">narg</span><span class="w"> </span><span class="o">&gt;</span><span class="w"> </span><span class="mi">8</span><span class="p">)</span><span class="w"> </span><span class="n">error</span><span class="o">-&gt;</span><span class="n">all</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;Incorrect args for pair coefficients&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="n">allocated</span><span class="p">)</span><span class="w"> </span><span class="n">allocate</span><span class="p">();</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">ilo</span><span class="p">,</span><span class="w"> </span><span class="n">ihi</span><span class="p">,</span><span class="w"> </span><span class="n">jlo</span><span class="p">,</span><span class="w"> </span><span class="n">jhi</span><span class="p">;</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">bounds</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">,</span><span class="w"> </span><span class="n">ilo</span><span class="p">,</span><span class="w"> </span><span class="n">ihi</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">bounds</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">,</span><span class="w"> </span><span class="n">jlo</span><span class="p">,</span><span class="w"> </span><span class="n">jhi</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">biga0_one</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">numeric</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span><span class="w"> </span><span class="nb">false</span><span class="p">,</span><span class="w"> </span><span class="n">lmp</span><span class="p">);</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">alpha_one</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">numeric</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="w"> </span><span class="nb">false</span><span class="p">,</span><span class="w"> </span><span class="n">lmp</span><span class="p">);</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">biga1_one</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">numeric</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="w"> </span><span class="nb">false</span><span class="p">,</span><span class="w"> </span><span class="n">lmp</span><span class="p">);</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">beta_one</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">numeric</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="w"> </span><span class="nb">false</span><span class="p">,</span><span class="w"> </span><span class="n">lmp</span><span class="p">);</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">r0_one</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">numeric</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">6</span><span class="p">],</span><span class="w"> </span><span class="nb">false</span><span class="p">,</span><span class="w"> </span><span class="n">lmp</span><span class="p">);</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">cut_one</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">cut_global</span><span class="p">;</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">narg</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">10</span><span class="p">)</span><span class="w"> </span><span class="n">cut_one</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">numeric</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">[</span><span class="mi">7</span><span class="p">],</span><span class="w"> </span><span class="nb">false</span><span class="p">,</span><span class="w"> </span><span class="n">lmp</span><span class="p">);</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">count</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">ilo</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">ihi</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">MAX</span><span class="p">(</span><span class="n">jlo</span><span class="p">,</span><span class="w"> </span><span class="n">i</span><span class="p">);</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">jhi</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga0_one</span><span class="p">;</span>
<span class="w"> </span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">alpha_one</span><span class="p">;</span>
<span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga1_one</span><span class="p">;</span>
<span class="w"> </span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">beta_one</span><span class="p">;</span>
<span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">r0_one</span><span class="p">;</span>
<span class="w"> </span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">cut_one</span><span class="p">;</span>
<span class="w"> </span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="w"> </span><span class="n">count</span><span class="o">++</span><span class="p">;</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">count</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="n">error</span><span class="o">-&gt;</span><span class="n">all</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;Incorrect args for pair coefficients&quot;</span><span class="p">);</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="initialization">
<h2>Initialization<a class="headerlink" href="#initialization" title="Link to this heading"></a></h2>
<p>The <code class="docutils literal notranslate"><span class="pre">init_one()</span></code> function is called during the <a class="reference internal" href="Developer_flow.html"><span class="doc">“init” phase</span></a> of a simulation. This is where potential parameters
are checked for completeness, derived parameters computed (e.g. the
“offset” of the potential energy at the cutoff distance for use with the
<a class="reference internal" href="pair_modify.html"><span class="doc">pair_modify shift yes</span></a> command). If a pair style
supports generating “mixed” parameters (i.e. where both atoms of a pair
have a different atom type) using a “mixing rule” from the parameters of
the type with itself, this is the place to compute and store those mixed
values. The <em>born/gauss</em> pair style does not support mixing, so we only
check for completeness. Another purpose of the <code class="docutils literal notranslate"><span class="pre">init_one()</span></code> function
is to symmetrize the potential parameter arrays. The return value of
the function is the cutoff for the given pair of atom types. This
information is used by the neighbor list code to determine the largest
cutoff and then build the neighbor lists accordingly.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> init for one type pair i,j and corresponding j,i</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">double</span><span class="w"> </span><span class="nf">PairBornGauss::init_one</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">j</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="n">error</span><span class="o">-&gt;</span><span class="n">all</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;All pair coeffs are not set&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">offset_flag</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="n">offset</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span>
<span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">])</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span><span class="w"> </span><span class="k">else</span>
<span class="w"> </span><span class="n">offset</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mf">0.0</span><span class="p">;</span>
<span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="n">alpha</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="n">beta</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="n">offset</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">offset</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">];</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="computing-forces-from-the-neighbor-list-required">
<h2>Computing forces from the neighbor list (required)<a class="headerlink" href="#computing-forces-from-the-neighbor-list-required" title="Link to this heading"></a></h2>
<p>The <code class="docutils literal notranslate"><span class="pre">compute()</span></code> function is the “workhorse” of a pair style. This is
where we have the nested loops over all pairs of particles from the
neighbor list to compute forces and - if needed - energies and virials.</p>
<p>The first part is to define some variables for later use and store
cached copies of data or pointers that we need to access frequently. Also,
this is a good place to call <code class="docutils literal notranslate"><span class="pre">Pair::ev_init()</span></code>, which initializes
several flags derived from the <cite>eflag</cite> and <cite>vflag</cite> parameters signaling
whether the energy and virial need to be tallied and whether only globally
or also per-atom.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ---------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::compute</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">eflag</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">vflag</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="w"> </span><span class="n">j</span><span class="p">,</span><span class="w"> </span><span class="n">ii</span><span class="p">,</span><span class="w"> </span><span class="n">jj</span><span class="p">,</span><span class="w"> </span><span class="n">inum</span><span class="p">,</span><span class="w"> </span><span class="n">jnum</span><span class="p">,</span><span class="w"> </span><span class="n">itype</span><span class="p">,</span><span class="w"> </span><span class="n">jtype</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">xtmp</span><span class="p">,</span><span class="w"> </span><span class="n">ytmp</span><span class="p">,</span><span class="w"> </span><span class="n">ztmp</span><span class="p">,</span><span class="w"> </span><span class="n">delx</span><span class="p">,</span><span class="w"> </span><span class="n">dely</span><span class="p">,</span><span class="w"> </span><span class="n">delz</span><span class="p">,</span><span class="w"> </span><span class="n">evdwl</span><span class="p">,</span><span class="w"> </span><span class="n">fpair</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">rsq</span><span class="p">,</span><span class="w"> </span><span class="n">r</span><span class="p">,</span><span class="w"> </span><span class="n">dr</span><span class="p">,</span><span class="w"> </span><span class="n">aexp</span><span class="p">,</span><span class="w"> </span><span class="n">bexp</span><span class="p">,</span><span class="w"> </span><span class="n">factor_lj</span><span class="p">;</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="o">*</span><span class="n">ilist</span><span class="p">,</span><span class="w"> </span><span class="o">*</span><span class="n">jlist</span><span class="p">,</span><span class="w"> </span><span class="o">*</span><span class="n">numneigh</span><span class="p">,</span><span class="w"> </span><span class="o">**</span><span class="n">firstneigh</span><span class="p">;</span>
<span class="w"> </span><span class="n">evdwl</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mf">0.0</span><span class="p">;</span>
<span class="w"> </span><span class="n">ev_init</span><span class="p">(</span><span class="n">eflag</span><span class="p">,</span><span class="w"> </span><span class="n">vflag</span><span class="p">);</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">**</span><span class="n">x</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">x</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">**</span><span class="n">f</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">f</span><span class="p">;</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="o">*</span><span class="n">type</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">type</span><span class="p">;</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">nlocal</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">nlocal</span><span class="p">;</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">*</span><span class="n">special_lj</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">force</span><span class="o">-&gt;</span><span class="n">special_lj</span><span class="p">;</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">newton_pair</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">force</span><span class="o">-&gt;</span><span class="n">newton_pair</span><span class="p">;</span>
<span class="w"> </span><span class="n">inum</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">list</span><span class="o">-&gt;</span><span class="n">inum</span><span class="p">;</span>
<span class="w"> </span><span class="n">ilist</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">list</span><span class="o">-&gt;</span><span class="n">ilist</span><span class="p">;</span>
<span class="w"> </span><span class="n">numneigh</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">list</span><span class="o">-&gt;</span><span class="n">numneigh</span><span class="p">;</span>
<span class="w"> </span><span class="n">firstneigh</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">list</span><span class="o">-&gt;</span><span class="n">firstneigh</span><span class="p">;</span>
</pre></div>
</div>
<p>The outer loop (index <em>i</em>) is over local atoms of our sub-domain.
Typically, the value of <cite>inum</cite> (the number of neighbor lists) is the
same as the number of local atoms (= atoms <em>owned</em> by this sub-domain).
But when the pair style is used as a sub-style of a <a class="reference internal" href="pair_hybrid.html"><span class="doc">hybrid pair
style</span></a> or neighbor list entries are removed with
<a class="reference internal" href="neigh_modify.html"><span class="doc">neigh_modify exclude</span></a>, this number may be
smaller. The array <code class="docutils literal notranslate"><span class="pre">list-&gt;ilist</span></code> has the (local) indices of the atoms
for which neighbor lists have been created. Then <code class="docutils literal notranslate"><span class="pre">list-&gt;numneigh</span></code> is
an <cite>inum</cite> sized array with the number of entries of each list of
neighbors, and <code class="docutils literal notranslate"><span class="pre">list-&gt;firstneigh</span></code> is a list of pointers to those lists.</p>
<p>For efficiency reasons, cached copies of some properties of the outer
loop atoms are also initialized.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="c1">// loop over neighbors of my atoms</span>
<span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">ii</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">ii</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">inum</span><span class="p">;</span><span class="w"> </span><span class="n">ii</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">ilist</span><span class="p">[</span><span class="n">ii</span><span class="p">];</span>
<span class="w"> </span><span class="n">xtmp</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">0</span><span class="p">];</span>
<span class="w"> </span><span class="n">ytmp</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">1</span><span class="p">];</span>
<span class="w"> </span><span class="n">ztmp</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">2</span><span class="p">];</span>
<span class="w"> </span><span class="n">itype</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">type</span><span class="p">[</span><span class="n">i</span><span class="p">];</span>
<span class="w"> </span><span class="n">jlist</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">firstneigh</span><span class="p">[</span><span class="n">i</span><span class="p">];</span>
<span class="w"> </span><span class="n">jnum</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">numneigh</span><span class="p">[</span><span class="n">i</span><span class="p">];</span>
</pre></div>
</div>
<p>The inner loop (index <em>j</em>) processes the neighbor lists. The neighbor
list code encodes extra information using the upper 3 bits. The 2
highest bits encode whether a pair is a regular pair of neighbor (= 0)
or a pair of 1-2 (= 1), 1-3 (= 2), or 1-4 (= 3) <a class="reference internal" href="special_bonds.html"><span class="doc">“special” neighbor</span></a>. The next highest bit encodes whether the pair stores
data in a <code class="docutils literal notranslate"><span class="pre">fix</span> <span class="pre">neigh/history</span></code> instance (an undocumented internal fix
style). The <code class="docutils literal notranslate"><span class="pre">sbmask()</span></code> inline function extracts those bits and
converts them into a number. This number is used to look up the
corresponding scaling factor for the non-bonded interaction from the
<code class="docutils literal notranslate"><span class="pre">force-&gt;special_lj</span></code> array and stores it in the <cite>factor_lj</cite> variable.
Due to the additional bits, the value of <em>j</em> would be out of range when
accessing data from per-atom arrays, so we apply the NEIGHMASK constant
with a bit-wise and operation to mask them out. This step <em>must</em> be
done, even if a pair style does not use special bond scaling of forces
and energies to avoid segmentation faults.</p>
<p>With the corrected <em>j</em> index, it is now possible to compute the distance
of the pair. For efficiency reasons, the square root is only taken
<em>after</em> the check for the cutoff (which has been stored as squared
cutoff by the <code class="docutils literal notranslate"><span class="pre">Pair</span></code> base class). For some pair styles, like the 12-6
Lennard-Jones potential, computing the square root can be avoided
entirely.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">jj</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">jj</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">jnum</span><span class="p">;</span><span class="w"> </span><span class="n">jj</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">jlist</span><span class="p">[</span><span class="n">jj</span><span class="p">];</span>
<span class="w"> </span><span class="n">factor_lj</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">special_lj</span><span class="p">[</span><span class="n">sbmask</span><span class="p">(</span><span class="n">j</span><span class="p">)];</span>
<span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&amp;=</span><span class="w"> </span><span class="n">NEIGHMASK</span><span class="p">;</span>
<span class="w"> </span><span class="n">delx</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">xtmp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">0</span><span class="p">];</span>
<span class="w"> </span><span class="n">dely</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">ytmp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">1</span><span class="p">];</span>
<span class="w"> </span><span class="n">delz</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">ztmp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">2</span><span class="p">];</span>
<span class="w"> </span><span class="n">rsq</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">delx</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">delx</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">dely</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dely</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">delz</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">delz</span><span class="p">;</span>
<span class="w"> </span><span class="n">jtype</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">type</span><span class="p">[</span><span class="n">j</span><span class="p">];</span>
</pre></div>
</div>
<p>The following block of code is the actual application of the model
potential to compute the force. Note, that <em>fpair</em> is the pair-wise
force divided by the distance, as this simplifies the projection of the
x-, y-, and z-components of the force vector by simply multiplying with
the respective distances in those directions.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">rsq</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">cutsq</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">])</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">r</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">sqrt</span><span class="p">(</span><span class="n">rsq</span><span class="p">);</span>
<span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">r</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">];</span>
<span class="w"> </span><span class="n">aexp</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">alpha</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">r</span><span class="p">);</span>
<span class="w"> </span><span class="n">bexp</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">beta</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="p">);</span>
<span class="w"> </span><span class="n">fpair</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">alpha</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">aexp</span><span class="p">;</span>
<span class="w"> </span><span class="n">fpair</span><span class="w"> </span><span class="o">-=</span><span class="w"> </span><span class="mf">2.0</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">beta</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">bexp</span><span class="p">;</span>
<span class="w"> </span><span class="n">fpair</span><span class="w"> </span><span class="o">*=</span><span class="w"> </span><span class="n">factor_lj</span><span class="w"> </span><span class="o">/</span><span class="w"> </span><span class="n">r</span><span class="p">;</span>
</pre></div>
</div>
<p>In the next block, the force is added to the per-atom force arrays. This
pair style uses a “half” neighbor list (each pair is listed only once)
so we take advantage of the fact that <span class="math notranslate nohighlight">\(\vec{F}_{ij} =
-\vec{F}_{ji}\)</span>, i.e. apply Newtons third law. The force is <em>always</em>
stored when the atom is a “local” atom. Index <em>i</em> atoms are always “local”
(i.e. <em>i</em> &lt; nlocal); index <em>j</em> atoms may be “ghost” atoms (<em>j</em> &gt;= nlocal).</p>
<p>Depending on the settings used with the <a class="reference internal" href="newton.html"><span class="doc">newton command</span></a>,
those pairs are only listed once globally (newton_pair == 1), then
forces must be stored even with ghost atoms and after all forces are
computed a “reverse communication” is performed to add those ghost atom
forces to their corresponding local atoms. If the setting is disabled,
then the extra communication is skipped, since for pairs straddling
sub-domain boundaries, the forces are computed twice and only stored
with the local atoms in the domain that <em>owns</em> it.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="n">f</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span><span class="w"> </span><span class="o">+=</span><span class="w"> </span><span class="n">delx</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">fpair</span><span class="p">;</span>
<span class="n">f</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span><span class="w"> </span><span class="o">+=</span><span class="w"> </span><span class="n">dely</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">fpair</span><span class="p">;</span>
<span class="n">f</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">2</span><span class="p">]</span><span class="w"> </span><span class="o">+=</span><span class="w"> </span><span class="n">delz</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">fpair</span><span class="p">;</span>
<span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">newton_pair</span><span class="w"> </span><span class="o">||</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">nlocal</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">f</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span><span class="w"> </span><span class="o">-=</span><span class="w"> </span><span class="n">delx</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">fpair</span><span class="p">;</span>
<span class="w"> </span><span class="n">f</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span><span class="w"> </span><span class="o">-=</span><span class="w"> </span><span class="n">dely</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">fpair</span><span class="p">;</span>
<span class="w"> </span><span class="n">f</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">2</span><span class="p">]</span><span class="w"> </span><span class="o">-=</span><span class="w"> </span><span class="n">delz</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">fpair</span><span class="p">;</span>
<span class="p">}</span>
</pre></div>
</div>
<p>The <code class="docutils literal notranslate"><span class="pre">ev_tally()</span></code> function tallies global or per-atom energy and
virial. For typical MD simulations, the potential energy is merely a
diagnostic and only needed on output. Similarly, the pressure may only
be computed for (infrequent) thermodynamic output. For all timesteps
where this information is not needed either, <cite>eflag</cite> or <cite>evflag</cite> are
zero and the computation and call to the tally function skipped. Note
that evdwl is initialized to zero at the beginning of the function, so
that it still is valid to access it, even if the energy is not computed
(e.g. when only the virial is needed).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">eflag</span><span class="p">)</span><span class="w"> </span><span class="n">evdwl</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">factor_lj</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="p">(</span><span class="n">aexp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">bexp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">offset</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]);</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">evflag</span><span class="p">)</span><span class="w"> </span><span class="n">ev_tally</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="w"> </span><span class="n">j</span><span class="p">,</span><span class="w"> </span><span class="n">nlocal</span><span class="p">,</span><span class="w"> </span><span class="n">newton_pair</span><span class="p">,</span><span class="w"> </span><span class="n">evdwl</span><span class="p">,</span><span class="w"> </span><span class="mf">0.0</span><span class="p">,</span><span class="w"> </span><span class="n">fpair</span><span class="p">,</span><span class="w"> </span><span class="n">delx</span><span class="p">,</span><span class="w"> </span><span class="n">dely</span><span class="p">,</span><span class="w"> </span><span class="n">delz</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">}</span>
</pre></div>
</div>
<p>If only the global virial is needed and no energy, then calls to
<code class="docutils literal notranslate"><span class="pre">ev_tally()</span></code> can be avoided altogether, and the global virial can be
computed more efficiently from the dot product of the total per-atom
force vector and the position vector of the corresponding atom,
<span class="math notranslate nohighlight">\(\vec{F}\cdot\vec{r}\)</span>. This has to be done <em>after</em> all pair-wise
forces are computed and <em>before</em> the reverse communication to collect
data from ghost atoms, since the position has to be the position that was
used to compute the force, i.e. <em>not</em> the “local” position if that ghost
atom is a periodic copy.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">vflag_fdotr</span><span class="p">)</span><span class="w"> </span><span class="n">virial_fdotr_compute</span><span class="p">();</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="computing-force-and-energy-for-a-single-pair">
<h2>Computing force and energy for a single pair<a class="headerlink" href="#computing-force-and-energy-for-a-single-pair" title="Link to this heading"></a></h2>
<p>Certain features in LAMMPS only require computing interactions between
individual pairs of atoms and the (optional) <code class="docutils literal notranslate"><span class="pre">single()</span></code> function is
needed to support those features (e.g. for tabulation of force and
energy with <a class="reference internal" href="pair_write.html"><span class="doc">pair_write</span></a>). This is a repetition of
the force kernel in the <code class="docutils literal notranslate"><span class="pre">compute()</span></code> function, but only for a single
pair of atoms, where the (squared) distance is provided as a parameter
(so it may not even be an existing distance between two specific atoms).
The energy is returned as the return value of the function and the force
as the <cite>fforce</cite> reference. Note, that this is, similar to how <em>fpair</em>
is used in the <code class="docutils literal notranslate"><span class="pre">compute()</span></code> function, the magnitude of the force along
the vector between the two atoms <em>divided</em> by the distance.</p>
<p>The <code class="docutils literal notranslate"><span class="pre">single()</span></code> function is optional, but it is expected to be
implemented for any true pair-wise additive potential. Many-body
potentials and special case potentials do not implement it. In a few
special cases (EAM, long-range Coulomb), the <code class="docutils literal notranslate"><span class="pre">single()</span></code> function
implements the pairwise additive part of the complete force interaction
and depends on either pre-computed properties (derivative of embedding
term for EAM) or post-computed non-pair-wise force contributions (KSpace
style in case of long-range Coulomb).</p>
<p>The member variable <cite>single_enable</cite> should be set to 0 in the
constructor, if it is not implemented (its default value is 1).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ---------------------------------------------------------------------- */</span>
<span class="kt">double</span><span class="w"> </span><span class="nf">PairBornGauss::single</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="cm">/*i*/</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="cm">/*j*/</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">itype</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">jtype</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">rsq</span><span class="p">,</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="cm">/*factor_coul*/</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">factor_lj</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">&amp;</span><span class="n">fforce</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">r</span><span class="p">,</span><span class="w"> </span><span class="n">dr</span><span class="p">,</span><span class="w"> </span><span class="n">aexp</span><span class="p">,</span><span class="w"> </span><span class="n">bexp</span><span class="p">;</span>
<span class="w"> </span><span class="n">r</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">sqrt</span><span class="p">(</span><span class="n">rsq</span><span class="p">);</span>
<span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">r</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">];</span>
<span class="w"> </span><span class="n">aexp</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">alpha</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">r</span><span class="p">);</span>
<span class="w"> </span><span class="n">bexp</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">beta</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="p">);</span>
<span class="w"> </span><span class="n">fforce</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">factor_lj</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="p">(</span><span class="n">alpha</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">aexp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="mf">2.0</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">dr</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">beta</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">bexp</span><span class="p">)</span><span class="w"> </span><span class="o">/</span><span class="w"> </span><span class="n">r</span><span class="p">;</span>
<span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="n">factor_lj</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="p">(</span><span class="n">aexp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">bexp</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">offset</span><span class="p">[</span><span class="n">itype</span><span class="p">][</span><span class="n">jtype</span><span class="p">]);</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="reading-and-writing-of-restart-files">
<h2>Reading and writing of restart files<a class="headerlink" href="#reading-and-writing-of-restart-files" title="Link to this heading"></a></h2>
<p>Support for writing and reading binary restart files is provided by the
following four functions. Writing is only done by MPI processor rank 0.
The output of global (not related to atom types) settings is usually
delegated to the <code class="docutils literal notranslate"><span class="pre">write_restart_settings()</span></code> function. This restart
facility is commonly only used, if there are small number of per-type
parameters. For potentials that use per-element parameters or tabulated
data and read these from files, those parameters and the name of the
potential file are not written to restart files and the <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff
command</span></a> has to re-issued when restarting. For pair styles
like “born/gauss” that do support writing to restart files, this is not
required.</p>
<p>Implementing the functions to read and write binary restart files is
optional. The member variable <cite>restartinfo</cite> should be set to 0 in the
constructor, if they are not implemented (its default value is 1).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> proc 0 writes to restart file</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::write_restart</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="n">fp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">write_restart_settings</span><span class="p">(</span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="w"> </span><span class="n">j</span><span class="p">;</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">i</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">int</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">])</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">}</span>
<span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> proc 0 writes to restart file</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::write_restart_settings</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="n">fp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">cut_global</span><span class="p">,</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">offset_flag</span><span class="p">,</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">int</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">fwrite</span><span class="p">(</span><span class="o">&amp;</span><span class="n">mix_flag</span><span class="p">,</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">int</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">);</span>
<span class="p">}</span>
</pre></div>
</div>
<p>Similarly, on reading, only MPI processor rank 0 has opened the restart
file and will read the data. The data is then distributed across all
parallel processes using calls to <code class="docutils literal notranslate"><span class="pre">MPI_Bcast()</span></code>. Before reading atom
type specific data, the corresponding storage needs to be allocated.
Order and number or storage size of items read must be exactly the same
as when writing, or else the data will be read incorrectly.</p>
<p>Reading uses the <a class="reference internal" href="Developer_utils.html#_CPPv4N9LAMMPS_NS5utils6sfreadEPKciPv6size_t6size_tP4FILEPKcP5Error" title="LAMMPS_NS::utils::sfread"><code class="xref cpp cpp-func docutils literal notranslate"><span class="pre">utils::sfread</span></code></a>
utility function to detect read errors and short reads, so that LAMMPS
can abort if that happens, e.g. when the restart file is corrupted.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> proc 0 reads from restart file, bcasts</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::read_restart</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="n">fp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">read_restart_settings</span><span class="p">(</span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">allocate</span><span class="p">();</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="w"> </span><span class="n">j</span><span class="p">;</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">me</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">comm</span><span class="o">-&gt;</span><span class="n">me</span><span class="p">;</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">i</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">me</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">int</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_INT</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">setflag</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">])</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">me</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_DOUBLE</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_DOUBLE</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_DOUBLE</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_DOUBLE</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_DOUBLE</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_DOUBLE</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">}</span>
<span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> proc 0 reads from restart file, bcasts</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::read_restart_settings</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="n">fp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">comm</span><span class="o">-&gt;</span><span class="n">me</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">cut_global</span><span class="p">,</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">double</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">offset_flag</span><span class="p">,</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">int</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="n">utils</span><span class="o">::</span><span class="n">sfread</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="w"> </span><span class="o">&amp;</span><span class="n">mix_flag</span><span class="p">,</span><span class="w"> </span><span class="k">sizeof</span><span class="p">(</span><span class="kt">int</span><span class="p">),</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="k">nullptr</span><span class="p">,</span><span class="w"> </span><span class="n">error</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">cut_global</span><span class="p">,</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_DOUBLE</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">offset_flag</span><span class="p">,</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_INT</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="w"> </span><span class="n">MPI_Bcast</span><span class="p">(</span><span class="o">&amp;</span><span class="n">mix_flag</span><span class="p">,</span><span class="w"> </span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">MPI_INT</span><span class="p">,</span><span class="w"> </span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">world</span><span class="p">);</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="writing-coefficients-to-data-files">
<h2>Writing coefficients to data files<a class="headerlink" href="#writing-coefficients-to-data-files" title="Link to this heading"></a></h2>
<p>The <code class="docutils literal notranslate"><span class="pre">write_data()</span></code> and <code class="docutils literal notranslate"><span class="pre">write_data_all()</span></code> functions are optional and
write out the current state of the <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff
settings</span></a> as “Pair Coeffs” or “PairIJ Coeffs” sections to a
data file when using the <a class="reference internal" href="write_data.html"><span class="doc">write_data command</span></a>. The
<code class="docutils literal notranslate"><span class="pre">write_data()</span></code> only writes out the diagonal elements of the pair
coefficient matrix, as that is required for the format of the “Pair
Coeffs” section. It is called when the “pair” option of the
<a class="reference internal" href="write_data.html"><span class="doc">write_data command</span></a> is “ii” (the default). This is
suitable for force fields where <em>all</em> off-diagonal terms of the pair
coefficient matrix are generated from mixing. If explicit settings for
off-diagonal elements were made, LAMMPS will print a warning, as those
would be lost. To avoid this, the “pair ij” option of <a class="reference internal" href="write_data.html"><span class="doc">write_data</span></a> can be used which will trigger calling the
<code class="docutils literal notranslate"><span class="pre">write_data_all()</span></code> function instead, which will write out all settings
of the pair coefficient matrix (regardless of whether they were
originally created from mixing or not).</p>
<p>These data file output functions are only useful for true pair-wise
additive potentials, where the potential parameters can be entered
through <em>multiple</em> <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff commands</span></a>. Pair styles
that require a single “pair_coeff * *” command are not compatible
with reading their parameters from data files. For pair styles like
<em>born/gauss</em> that do support writing to data files, the potential
parameters will be read from the data file, if present, and
<a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff commands</span></a> may not be needed.</p>
<p>The member variable <code class="docutils literal notranslate"><span class="pre">writedata</span></code> should be set to 1 in the constructor,
if these functions are implemented (the default value is 0).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> proc 0 writes to data file</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::write_data</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="n">fp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span>
<span class="w"> </span><span class="n">fprintf</span><span class="p">(</span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;%d %g %g %g %g %g</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">,</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">i</span><span class="p">],</span><span class="w"> </span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">i</span><span class="p">],</span><span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">i</span><span class="p">],</span><span class="w"> </span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">i</span><span class="p">],</span>
<span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">i</span><span class="p">]);</span>
<span class="p">}</span>
<span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> proc 0 writes all pairs to data file</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairBornGauss::write_data_all</span><span class="p">(</span><span class="kt">FILE</span><span class="w"> </span><span class="o">*</span><span class="n">fp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">i</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">&lt;=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">ntypes</span><span class="p">;</span><span class="w"> </span><span class="n">j</span><span class="o">++</span><span class="p">)</span>
<span class="w"> </span><span class="n">fprintf</span><span class="p">(</span><span class="n">fp</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;%d %d %g %g %g %g %g %g</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">,</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="w"> </span><span class="n">j</span><span class="p">,</span><span class="w"> </span><span class="n">biga0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="n">alpha</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="n">biga1</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span>
<span class="w"> </span><span class="n">beta</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="n">r0</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">],</span><span class="w"> </span><span class="n">cut</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]);</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="give-access-to-internal-data">
<h2>Give access to internal data<a class="headerlink" href="#give-access-to-internal-data" title="Link to this heading"></a></h2>
<p>The purpose of the <code class="docutils literal notranslate"><span class="pre">extract()</span></code> function is to facilitate access to
internal data of the pair style by other parts of LAMMPS. One possible
application is to use <a class="reference internal" href="fix_adapt.html"><span class="doc">fix adapt</span></a> to gradually change
potential parameters during a run. Here, we implement access to the
pair coefficient matrix parameters.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ---------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="o">*</span><span class="nf">PairBornGauss::extract</span><span class="p">(</span><span class="k">const</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">*</span><span class="n">str</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="o">&amp;</span><span class="n">dim</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">dim</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">2</span><span class="p">;</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">strcmp</span><span class="p">(</span><span class="n">str</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;biga0&quot;</span><span class="p">)</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="p">(</span><span class="kt">void</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="n">biga0</span><span class="p">;</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">strcmp</span><span class="p">(</span><span class="n">str</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;biga1&quot;</span><span class="p">)</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="p">(</span><span class="kt">void</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="n">biga1</span><span class="p">;</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">strcmp</span><span class="p">(</span><span class="n">str</span><span class="p">,</span><span class="w"> </span><span class="s">&quot;r0&quot;</span><span class="p">)</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="p">(</span><span class="kt">void</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="n">r0</span><span class="p">;</span>
<span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="k">nullptr</span><span class="p">;</span>
<span class="p">}</span>
</pre></div>
</div>
<p>Since the mercury potential, for which we have implemented the
born/gauss pair style, has a temperature dependent parameter “biga1”, we
can automatically adapt the potential based on the Taylor-MacLaurin
expansion for “biga1” when performing a simulation with a temperature
ramp. LAMMPS commands for that application are given below:</p>
<div class="highlight-LAMMPS notranslate"><div class="highlight"><pre><span></span><span class="k">variable </span><span class="nv nv-Identifier">tlo</span><span class="w"> </span><span class="n">index</span><span class="w"> </span><span class="m">300.0</span>
<span class="k">variable </span><span class="nv nv-Identifier">thi</span><span class="w"> </span><span class="n">index</span><span class="w"> </span><span class="m">600.0</span>
<span class="k">variable </span><span class="nv nv-Identifier">temp</span><span class="w"> </span><span class="n">equal</span><span class="w"> </span><span class="n">ramp</span><span class="nv">(v_tlo,v_thi)</span>
<span class="k">variable </span><span class="nv nv-Identifier">biga1</span><span class="w"> </span><span class="n">equal</span><span class="w"> </span><span class="nv">(-2.58717e-8*v_temp+8.40841e-5)</span><span class="o">*</span><span class="n">v_temp</span><span class="o">+</span><span class="m">1.97475e-2</span>
<span class="k">fix </span><span class="nv nv-Identifier">1</span><span class="w"> </span><span class="nv nv-Identifier">all</span><span class="w"> </span><span class="n">nvt</span><span class="w"> </span><span class="n">temp</span><span class="w"> </span><span class="nv">${tlo}</span><span class="w"> </span><span class="nv">${thi}</span><span class="w"> </span><span class="m">0.1</span>
<span class="k">fix </span><span class="nv nv-Identifier">2</span><span class="w"> </span><span class="nv nv-Identifier">all</span><span class="w"> </span><span class="n">adapt</span><span class="w"> </span><span class="m">1</span><span class="w"> </span><span class="n">pair</span><span class="w"> </span><span class="n">born</span><span class="o">/</span><span class="n">gauss</span><span class="w"> </span><span class="n">biga1</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">v_biga1</span>
</pre></div>
</div>
</section>
</section>
<section id="case-2-a-many-body-potential">
<h1><span class="section-number">4.8.4. </span>Case 2: a many-body potential<a class="headerlink" href="#case-2-a-many-body-potential" title="Link to this heading"></a></h1>
<p>Since there is a detailed description of the purpose and general layout
of a pair style in the previous case, we will focus on where the
implementation of a typical many-body potential <em>differs</em> from a
pair-wise additive potential. We will use the implementation of the
Tersoff potential as <a class="reference internal" href="pair_tersoff.html"><span class="doc">pair_style tersoff</span></a> as an
example. The complete implementation can be found in the files
<code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_tersoff.cpp</span></code> and <code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_tersoff.h</span></code> of
the LAMMPS source code.</p>
<section id="constructor">
<h2>Constructor<a class="headerlink" href="#constructor" title="Link to this heading"></a></h2>
<p>In the constructor, several <a class="reference internal" href="Modify_pair.html"><span class="doc">pair style flags</span></a> must
be set differently for many-body potentials:</p>
<ul class="simple">
<li><p>the potential is not pair-wise additive, so the <code class="docutils literal notranslate"><span class="pre">single()</span></code> function
cannot be used. This is indicated by setting the <cite>single_enable</cite>
member variable to 0 (default value is 1)</p></li>
<li><p>many-body potentials are usually not written to <a class="reference internal" href="write_restart.html"><span class="doc">binary
restart files</span></a>. This is indicated by setting the member
variable <cite>restartinfo</cite> to 0 (default is 1)</p></li>
<li><p>many-body potentials typically read <em>all</em> parameters from a file which
stores parameters indexed with a string (e.g. the element). For this,
only a single <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff * *</span></a> command is allowed.
This requirement is set and checked for, when the member variable
<cite>one_coeff</cite> is set to 1 (default value is 0)</p></li>
<li><p>many-body potentials can produce incorrect results if pairs of atoms
are excluded from the neighbor list, e.g. explicitly by
<a class="reference internal" href="neigh_modify.html"><span class="doc">neigh_modify exclude</span></a> or implicitly through
defining bonds, angles, etc. and having a <a class="reference internal" href="special_bonds.html"><span class="doc">special_bonds setting</span></a> that is not “special_bonds lj/coul 1.0 1.0 1.0”.
LAMMPS will check for this and print a suitable warning, when the
member variable <cite>manybody_flag</cite> is set to 1 (default value is 0).</p></li>
</ul>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="n">PairTersoff</span><span class="o">::</span><span class="n">PairTersoff</span><span class="p">(</span><span class="n">LAMMPS</span><span class="w"> </span><span class="o">*</span><span class="n">lmp</span><span class="p">)</span><span class="w"> </span><span class="o">:</span><span class="w"> </span><span class="n">Pair</span><span class="p">(</span><span class="n">lmp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">single_enable</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="n">restartinfo</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="n">one_coeff</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="w"> </span><span class="n">manybody_flag</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
</pre></div>
</div>
</section>
<section id="neighbor-list-request">
<h2>Neighbor list request<a class="headerlink" href="#neighbor-list-request" title="Link to this heading"></a></h2>
<p>For computing the three-body interactions of the Tersoff potential a
“full” neighbor list (both atoms of a pair are listed in each others
neighbor list) is required. By default a “half” neighbor list is
requested (each pair is listed only once). The request is made in
the <code class="docutils literal notranslate"><span class="pre">init_style()</span></code> function. A more in-depth discussion of neighbor
lists in LAMMPS and how to request them is in <a class="reference internal" href="Developer_notes.html#request-neighbor-list"><span class="std std-ref">this section of the
documentation</span></a></p>
<p>Also, additional conditions must be met for some global settings which
are checked in the <code class="docutils literal notranslate"><span class="pre">init_style()</span></code> function.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="cm">/* ----------------------------------------------------------------------</span>
<span class="cm"> init specific to this pair style</span>
<span class="cm">------------------------------------------------------------------------- */</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairTersoff::init_style</span><span class="p">()</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">tag_enable</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span>
<span class="w"> </span><span class="n">error</span><span class="o">-&gt;</span><span class="n">all</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="s">&quot;Pair style Tersoff requires atom IDs&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">force</span><span class="o">-&gt;</span><span class="n">newton_pair</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span>
<span class="w"> </span><span class="n">error</span><span class="o">-&gt;</span><span class="n">all</span><span class="p">(</span><span class="n">FLERR</span><span class="p">,</span><span class="s">&quot;Pair style Tersoff requires newton pair on&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="c1">// need a full neighbor list</span>
<span class="w"> </span><span class="n">neighbor</span><span class="o">-&gt;</span><span class="n">add_request</span><span class="p">(</span><span class="k">this</span><span class="p">,</span><span class="n">NeighConst</span><span class="o">::</span><span class="n">REQ_FULL</span><span class="p">);</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="computing-forces-from-the-neighbor-list">
<h2>Computing forces from the neighbor list<a class="headerlink" href="#computing-forces-from-the-neighbor-list" title="Link to this heading"></a></h2>
<p>Computing forces for a many-body potential is usually more complex than
for a pair-wise additive potential and there are multiple components.
For Tersoff, there is a pair-wise additive two-body term (two nested
loops over indices <em>i</em> and <em>j</em>) and a three-body term (three nested
loops over indices <em>i</em>, <em>j</em>, and <em>k</em>). Since the neighbor list has
all neighbors up to the maximum cutoff (for the two-body term), but
the three-body interactions have a significantly shorter cutoff,
a “short neighbor list” is also constructed at the same time while computing
the two-body term and looping over the neighbor list for the first time.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">rsq</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">cutshortsq</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">neighshort</span><span class="p">[</span><span class="n">numshort</span><span class="o">++</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">j</span><span class="p">;</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">numshort</span><span class="w"> </span><span class="o">&gt;=</span><span class="w"> </span><span class="n">maxshort</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">maxshort</span><span class="w"> </span><span class="o">+=</span><span class="w"> </span><span class="n">maxshort</span><span class="o">/</span><span class="mi">2</span><span class="p">;</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">grow</span><span class="p">(</span><span class="n">neighshort</span><span class="p">,</span><span class="n">maxshort</span><span class="p">,</span><span class="s">&quot;pair:neighshort&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">}</span>
</pre></div>
</div>
<p>For the two-body term, only a half neighbor list would be needed, even
though we have requested a full list (for the three-body loops).
Rather than computing all interactions twice, we skip over half of
the entries. This is done in a slightly complex way to make certain
the same choice is made across all subdomains and so that there is
no load imbalance introduced.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="n">jtag</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">tag</span><span class="p">[</span><span class="n">j</span><span class="p">];</span>
<span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">itag</span><span class="w"> </span><span class="o">&gt;</span><span class="w"> </span><span class="n">jtag</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">((</span><span class="n">itag</span><span class="o">+</span><span class="n">jtag</span><span class="p">)</span><span class="w"> </span><span class="o">%</span><span class="w"> </span><span class="mi">2</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span><span class="w"> </span><span class="k">continue</span><span class="p">;</span>
<span class="p">}</span><span class="w"> </span><span class="k">else</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">itag</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">jtag</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">((</span><span class="n">itag</span><span class="o">+</span><span class="n">jtag</span><span class="p">)</span><span class="w"> </span><span class="o">%</span><span class="w"> </span><span class="mi">2</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">1</span><span class="p">)</span><span class="w"> </span><span class="k">continue</span><span class="p">;</span>
<span class="p">}</span><span class="w"> </span><span class="k">else</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">2</span><span class="p">]</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">2</span><span class="p">])</span><span class="w"> </span><span class="k">continue</span><span class="p">;</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">2</span><span class="p">]</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="n">ztmp</span><span class="w"> </span><span class="o">&amp;&amp;</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">ytmp</span><span class="p">)</span><span class="w"> </span><span class="k">continue</span><span class="p">;</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">2</span><span class="p">]</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="n">ztmp</span><span class="w"> </span><span class="o">&amp;&amp;</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="n">ytmp</span><span class="w"> </span><span class="o">&amp;&amp;</span><span class="w"> </span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">xtmp</span><span class="p">)</span><span class="w"> </span><span class="k">continue</span><span class="p">;</span>
<span class="p">}</span>
</pre></div>
</div>
<p>For the three-body term, there is one additional nested loop and it uses
the “short” neighbor list, accumulated previously.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="c1">// three-body interactions</span>
<span class="c1">// skip immediately if I-J is not within cutoff</span>
<span class="kt">double</span><span class="w"> </span><span class="n">fjxtmp</span><span class="p">,</span><span class="n">fjytmp</span><span class="p">,</span><span class="n">fjztmp</span><span class="p">;</span>
<span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">jj</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">jj</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">numshort</span><span class="p">;</span><span class="w"> </span><span class="n">jj</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">neighshort</span><span class="p">[</span><span class="n">jj</span><span class="p">];</span>
<span class="w"> </span><span class="n">jtype</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">map</span><span class="p">[</span><span class="n">type</span><span class="p">[</span><span class="n">j</span><span class="p">]];</span>
<span class="w"> </span><span class="p">[...]</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">kk</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">kk</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">numshort</span><span class="p">;</span><span class="w"> </span><span class="n">kk</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">jj</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="n">kk</span><span class="p">)</span><span class="w"> </span><span class="k">continue</span><span class="p">;</span>
<span class="w"> </span><span class="n">k</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">neighshort</span><span class="p">[</span><span class="n">kk</span><span class="p">];</span>
<span class="w"> </span><span class="n">ktype</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">map</span><span class="p">[</span><span class="n">type</span><span class="p">[</span><span class="n">k</span><span class="p">]];</span>
<span class="w"> </span><span class="p">[...]</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">[...]</span>
</pre></div>
</div>
</section>
<section id="reading-potential-parameters">
<h2>Reading potential parameters<a class="headerlink" href="#reading-potential-parameters" title="Link to this heading"></a></h2>
<p>For the Tersoff potential, the parameters are listed in a file and
associated with triples of elements. Because we have set the
<code class="docutils literal notranslate"><span class="pre">one_coeff</span></code> flag to 1 in the constructor, there may only be a single
<a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff * *</span></a> line in the input for this pair
style, and as a consequence the <code class="docutils literal notranslate"><span class="pre">coeff()</span></code> function will only be called
once. Thus, the <code class="docutils literal notranslate"><span class="pre">coeff()</span></code> function has to do three tasks, each of
which is delegated to a function in the <code class="docutils literal notranslate"><span class="pre">PairTersoff</span></code> class:</p>
<ol class="arabic simple">
<li><p>map elements to atom types. Those follow the potential file name in the
command arguments and are processed by the <code class="docutils literal notranslate"><span class="pre">map_element2type()</span></code> function.</p></li>
<li><p>read and parse the potential parameter file in the <code class="docutils literal notranslate"><span class="pre">read_file()</span></code> function.</p></li>
<li><p>Build data structures where the original and derived parameters are
indexed by all possible triples of atom types and thus can be looked
up quickly in the loops for the force computation</p></li>
</ol>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="kt">void</span><span class="w"> </span><span class="nf">PairTersoff::coeff</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">narg</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="n">arg</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="n">allocated</span><span class="p">)</span><span class="w"> </span><span class="n">allocate</span><span class="p">();</span>
<span class="w"> </span><span class="n">map_element2type</span><span class="p">(</span><span class="n">narg</span><span class="mi">-3</span><span class="p">,</span><span class="n">arg</span><span class="o">+</span><span class="mi">3</span><span class="p">);</span>
<span class="w"> </span><span class="c1">// read potential file and initialize potential parameters</span>
<span class="w"> </span><span class="n">read_file</span><span class="p">(</span><span class="n">arg</span><span class="p">[</span><span class="mi">2</span><span class="p">]);</span>
<span class="w"> </span><span class="n">setup_params</span><span class="p">();</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
</section>
<section id="case-3-a-potential-requiring-communication">
<h1><span class="section-number">4.8.5. </span>Case 3: a potential requiring communication<a class="headerlink" href="#case-3-a-potential-requiring-communication" title="Link to this heading"></a></h1>
<p>For some models, the interactions between atoms depends on properties of
their environment which have to be computed <em>before</em> the the forces can
be computed. Since LAMMPS is designed to run in parallel using a
<a class="reference internal" href="Developer_par_part.html"><span class="doc">domain decomposition strategy</span></a>, not all
information of the atoms may be directly available and thus
communication steps may be need to collect data from ghost atoms of
neighboring subdomains or send data to ghost atoms for application
during the pairwise computation.</p>
<p>Specifically, two communication patterns are needed: a “reverse
communication” and a “forward communication”. The reverse communication
collects data added to “ghost” atoms from neighboring sub-domains and
sums it to their corresponding “local” atoms. This communication is
only required and thus executed when the <code class="docutils literal notranslate"><span class="pre">Force::newton_pair</span></code> setting
is 1 (i.e. <a class="reference internal" href="newton.html"><span class="doc">newton on</span></a>, the default). The forward
communication is used to copy computed per-atom data from “local” atoms
to their corresponding “ghost” atoms in neighboring sub-domains.</p>
<p>For this we will look at how the embedding term of the <a class="reference internal" href="pair_eam.html"><span class="doc">embedded
atom potential EAM</span></a> is implemented in LAMMPS. The complete
implementation of this pair style can be found in the files
<code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_eam.cpp</span></code> and <code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_eam.h</span></code> of the
LAMMPS source code.</p>
<section id="allocating-additional-per-atom-storage">
<h2>Allocating additional per-atom storage<a class="headerlink" href="#allocating-additional-per-atom-storage" title="Link to this heading"></a></h2>
<p>First suitable (local) per-atom arrays (<cite>rho</cite>, <cite>fp</cite>, <cite>numforce</cite>) are
allocated. These have to be large enough to include ghost atoms, are not
used outside the <code class="docutils literal notranslate"><span class="pre">compute()</span></code> function and are re-initialized to zero
once per timestep.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">nmax</span><span class="w"> </span><span class="o">&gt;</span><span class="w"> </span><span class="n">nmax</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">rho</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">fp</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">destroy</span><span class="p">(</span><span class="n">numforce</span><span class="p">);</span>
<span class="w"> </span><span class="n">nmax</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">atom</span><span class="o">-&gt;</span><span class="n">nmax</span><span class="p">;</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">rho</span><span class="p">,</span><span class="n">nmax</span><span class="p">,</span><span class="s">&quot;pair:rho&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">fp</span><span class="p">,</span><span class="n">nmax</span><span class="p">,</span><span class="s">&quot;pair:fp&quot;</span><span class="p">);</span>
<span class="w"> </span><span class="n">memory</span><span class="o">-&gt;</span><span class="n">create</span><span class="p">(</span><span class="n">numforce</span><span class="p">,</span><span class="n">nmax</span><span class="p">,</span><span class="s">&quot;pair:numforce&quot;</span><span class="p">);</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="reverse-communication">
<h2>Reverse communication<a class="headerlink" href="#reverse-communication" title="Link to this heading"></a></h2>
<p>Then a first loop over all pairs (<em>i</em> and <em>j</em>) is performed, where data
is stored in the <cite>rho</cite> array representing the electron density at the site of
<em>i</em> contributed from all neighbors <em>j</em>. Since the EAM pair style uses
a half neighbor list (for efficiency reasons), a reverse communication is
needed to collect the contributions to <cite>rho</cite> from ghost atoms (only if
<a class="reference internal" href="newton.html"><span class="doc">newton on</span></a> is set for pair styles).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">newton_pair</span><span class="p">)</span><span class="w"> </span><span class="n">comm</span><span class="o">-&gt;</span><span class="n">reverse_comm</span><span class="p">(</span><span class="k">this</span><span class="p">);</span>
</pre></div>
</div>
<p>To support the reverse communication, two functions must be defined:
<code class="docutils literal notranslate"><span class="pre">pack_reverse_comm()</span></code> that copies relevant data into a buffer for ghost
atoms and <code class="docutils literal notranslate"><span class="pre">unpack_reverse_comm()</span></code> that takes the collected data and adds
it to the <cite>rho</cite> array for the corresponding local atoms that match the
ghost atoms. In order to allocate sufficiently sized buffers, a flag
must be set in the pair style constructor. Since in this case a single
double precision number is communicated per atom, the <cite>comm_reverse</cite>
member variable is set to 1 (default is 0 = no reverse communication).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="kt">int</span><span class="w"> </span><span class="nf">PairEAM::pack_reverse_comm</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">first</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">*</span><span class="n">buf</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="n">m</span><span class="p">,</span><span class="n">last</span><span class="p">;</span>
<span class="w"> </span><span class="n">m</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="n">last</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">first</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">n</span><span class="p">;</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">first</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">last</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="n">buf</span><span class="p">[</span><span class="n">m</span><span class="o">++</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">rho</span><span class="p">[</span><span class="n">i</span><span class="p">];</span>
<span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="n">m</span><span class="p">;</span>
<span class="p">}</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairEAM::unpack_reverse_comm</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="o">*</span><span class="n">list</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">*</span><span class="n">buf</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">m</span><span class="p">;</span>
<span class="w"> </span><span class="n">m</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">n</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">list</span><span class="p">[</span><span class="n">i</span><span class="p">];</span>
<span class="w"> </span><span class="n">rho</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">+=</span><span class="w"> </span><span class="n">buf</span><span class="p">[</span><span class="n">m</span><span class="o">++</span><span class="p">];</span>
<span class="w"> </span><span class="p">}</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
<section id="forward-communication">
<h2>Forward communication<a class="headerlink" href="#forward-communication" title="Link to this heading"></a></h2>
<p>From the density array <cite>rho</cite>, the derivative of the embedding energy
<cite>fp</cite> is computed. The computation is only done for “local” atoms, but
for the force computation, that property also is needed on ghost atoms.
For that a forward communication is needed.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="n">comm</span><span class="o">-&gt;</span><span class="n">forward_comm</span><span class="p">(</span><span class="k">this</span><span class="p">);</span>
</pre></div>
</div>
<p>Similar to the reverse communication, this requires implementing a
<code class="docutils literal notranslate"><span class="pre">pack_forward_comm()</span></code> and an <code class="docutils literal notranslate"><span class="pre">unpack_forward_comm()</span></code> function.
Since there is one double precision number per atom that needs to be
communicated, we must set the <cite>comm_forward</cite> member variable to 1
(default is 0 = no forward communication).</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="kt">int</span><span class="w"> </span><span class="nf">PairEAM::pack_forward_comm</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="o">*</span><span class="n">list</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">*</span><span class="n">buf</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">pbc_flag</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="o">*</span><span class="n">pbc</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">m</span><span class="p">;</span>
<span class="w"> </span><span class="n">m</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">n</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="n">j</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">list</span><span class="p">[</span><span class="n">i</span><span class="p">];</span>
<span class="w"> </span><span class="n">buf</span><span class="p">[</span><span class="n">m</span><span class="o">++</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">fp</span><span class="p">[</span><span class="n">j</span><span class="p">];</span>
<span class="w"> </span><span class="p">}</span>
<span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="n">m</span><span class="p">;</span>
<span class="p">}</span>
<span class="kt">void</span><span class="w"> </span><span class="nf">PairEAM::unpack_forward_comm</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">first</span><span class="p">,</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="o">*</span><span class="n">buf</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="kt">int</span><span class="w"> </span><span class="n">i</span><span class="p">,</span><span class="n">m</span><span class="p">,</span><span class="n">last</span><span class="p">;</span>
<span class="w"> </span><span class="n">m</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span>
<span class="w"> </span><span class="n">last</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">first</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">n</span><span class="p">;</span>
<span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">first</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">&lt;</span><span class="w"> </span><span class="n">last</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="o">++</span><span class="p">)</span><span class="w"> </span><span class="n">fp</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">buf</span><span class="p">[</span><span class="n">m</span><span class="o">++</span><span class="p">];</span>
<span class="p">}</span>
</pre></div>
</div>
</section>
</section>
<section id="case-4-potentials-without-a-compute-function">
<h1><span class="section-number">4.8.6. </span>Case 4: potentials without a compute() function<a class="headerlink" href="#case-4-potentials-without-a-compute-function" title="Link to this heading"></a></h1>
<p>A small number of pair style classes do not implement a <code class="docutils literal notranslate"><span class="pre">compute()</span></code>
function, but instead use that of a different pair style.</p>
<section id="embedded-atom-variants-eam-fs-and-eam-alloy">
<h2>Embedded atom variants “eam/fs” and “eam/alloy”<a class="headerlink" href="#embedded-atom-variants-eam-fs-and-eam-alloy" title="Link to this heading"></a></h2>
<p>The pair styles <a class="reference internal" href="pair_eam.html"><span class="doc">eam/fs and eam/alloy</span></a> share the same
model and potential function as the <a class="reference internal" href="pair_eam.html"><span class="doc">eam pair style</span></a>.
They differ in the format of the potential files. Pair style <a class="reference internal" href="pair_eam.html"><span class="doc">eam</span></a> supports only potential files for single elements. For
multi-element systems, the mixed terms are computed from mixed
parameters. The <em>eam/fs</em> and <em>eam/alloy</em> pair styles, however,
<strong>require</strong> the use of a single potential file for all elements where
the mixed element potential is included in the tabulation. That enables
more accurate models for alloys, since the mixed terms can be adjusted
for a better representation of material properties compared to terms
created from mixing of per-element terms in the <code class="docutils literal notranslate"><span class="pre">PairEAM</span></code> class.</p>
<p>We take a closer at the <em>eam/alloy</em> pair style. The complete
implementation is in the files <code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_eam_alloy.cpp</span></code> and
<code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_eam_alloy.h</span></code>.</p>
<p>The <code class="docutils literal notranslate"><span class="pre">PairEAMAlloy</span></code> class is derived from <code class="docutils literal notranslate"><span class="pre">PairEAM</span></code> and not <code class="docutils literal notranslate"><span class="pre">Pair</span></code>
and overrides only a small number of functions:</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="k">class</span><span class="w"> </span><span class="nc">PairEAMAlloy</span><span class="w"> </span><span class="o">:</span><span class="w"> </span><span class="k">virtual</span><span class="w"> </span><span class="k">public</span><span class="w"> </span><span class="n">PairEAM</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">public</span><span class="o">:</span>
<span class="w"> </span><span class="n">PairEAMAlloy</span><span class="p">(</span><span class="k">class</span><span class="w"> </span><span class="nc">LAMMPS</span><span class="w"> </span><span class="o">*</span><span class="p">);</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">coeff</span><span class="p">(</span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="k">protected</span><span class="o">:</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="n">read_file</span><span class="p">(</span><span class="kt">char</span><span class="w"> </span><span class="o">*</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">file2array</span><span class="p">()</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="p">};</span>
</pre></div>
</div>
<p>All other functionality is inherited from the base classes. In the
constructor we set the <code class="docutils literal notranslate"><span class="pre">one_coeff</span></code> flag and the <code class="docutils literal notranslate"><span class="pre">many_body</span></code> flag to
1 to indicate the different behavior.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="n">PairEAMAlloy</span><span class="o">::</span><span class="n">PairEAMAlloy</span><span class="p">(</span><span class="n">LAMMPS</span><span class="w"> </span><span class="o">*</span><span class="n">lmp</span><span class="p">)</span><span class="w"> </span><span class="o">:</span><span class="w"> </span><span class="n">PairEAM</span><span class="p">(</span><span class="n">lmp</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">one_coeff</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="w"> </span><span class="n">manybody_flag</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="p">}</span>
</pre></div>
</div>
<p>The <code class="docutils literal notranslate"><span class="pre">coeff()</span></code> function (not shown here) implements the different
behavior when processing the <a class="reference internal" href="pair_coeff.html"><span class="doc">pair_coeff command</span></a>.
The <code class="docutils literal notranslate"><span class="pre">read_file()</span></code> and <code class="docutils literal notranslate"><span class="pre">file2array()</span></code> replace the corresponding
<code class="docutils literal notranslate"><span class="pre">PairEAM</span></code> class functions to accommodate the different data and
file format.</p>
</section>
<section id="airebo-and-airebo-m-potentials">
<h2>AIREBO and AIREBO-M potentials<a class="headerlink" href="#airebo-and-airebo-m-potentials" title="Link to this heading"></a></h2>
<p>The AIREBO-M potential differs from the better known AIREBO potential in
that it use a Morse potential instead of a Lennard-Jones potential for
non-bonded interactions. Since this difference is very minimal compared
to the entire potential, both potentials are implemented in the
<code class="docutils literal notranslate"><span class="pre">PairAIREBO</span></code> class and which non-bonded potential is used is
determined by the value of the <code class="docutils literal notranslate"><span class="pre">morseflag</span></code> flag, which would be set to
either 0 or 1.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="k">class</span><span class="w"> </span><span class="nc">PairAIREBOMorse</span><span class="w"> </span><span class="o">:</span><span class="w"> </span><span class="k">public</span><span class="w"> </span><span class="n">PairAIREBO</span><span class="w"> </span><span class="p">{</span>
<span class="w"> </span><span class="k">public</span><span class="o">:</span>
<span class="w"> </span><span class="n">PairAIREBOMorse</span><span class="p">(</span><span class="k">class</span><span class="w"> </span><span class="nc">LAMMPS</span><span class="w"> </span><span class="o">*</span><span class="p">);</span>
<span class="w"> </span><span class="kt">void</span><span class="w"> </span><span class="nf">settings</span><span class="p">(</span><span class="kt">int</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="p">)</span><span class="w"> </span><span class="k">override</span><span class="p">;</span>
<span class="p">};</span>
</pre></div>
</div>
<p>The <code class="docutils literal notranslate"><span class="pre">morseflag</span></code> variable defaults to 0 and is set to 1 in the
<code class="docutils literal notranslate"><span class="pre">PairAIREBOMorse::settings()</span></code> function which is called by the
<a class="reference internal" href="pair_style.html"><span class="doc">pair_style</span></a> command. This function delegates all
command argument processing and setting of other parameters to the
<code class="docutils literal notranslate"><span class="pre">PairAIREBO::settings()</span></code> function of the base class.</p>
<div class="highlight-c++ notranslate"><div class="highlight"><pre><span></span><span class="kt">void</span><span class="w"> </span><span class="nf">PairAIREBOMorse::settings</span><span class="p">(</span><span class="kt">int</span><span class="w"> </span><span class="n">narg</span><span class="p">,</span><span class="w"> </span><span class="kt">char</span><span class="w"> </span><span class="o">**</span><span class="n">arg</span><span class="p">)</span>
<span class="p">{</span>
<span class="w"> </span><span class="n">PairAIREBO</span><span class="o">::</span><span class="n">settings</span><span class="p">(</span><span class="n">narg</span><span class="p">,</span><span class="w"> </span><span class="n">arg</span><span class="p">);</span>
<span class="w"> </span><span class="n">morseflag</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="p">}</span>
</pre></div>
</div>
<p>The complete implementation is in the files
<code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_airebo.cpp</span></code>, <code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_airebo.h</span></code>,
<code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_airebo_morse.cpp</span></code>,
<code class="docutils literal notranslate"><span class="pre">src/MANYBODY/pair_airebo_morse.h</span></code>.</p>
<hr class="docutils" />
<p id="bomont2"><strong>(Bomont)</strong> Bomont, Bretonnet, J. Chem. Phys. 124, 054504 (2006)</p>
</section>
</section>
</div>
</div>
<footer><div class="rst-footer-buttons" role="navigation" aria-label="Footer">
<a href="Developer_write.html" class="btn btn-neutral float-left" title="4.8. Writing new styles" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left" aria-hidden="true"></span> Previous</a>
<a href="Developer_write_fix.html" class="btn btn-neutral float-right" title="4.8.7. Writing a new fix style" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right" aria-hidden="true"></span></a>
</div>
<hr/>
<div role="contentinfo">
<p>&#169; Copyright 2003-2025 Sandia Corporation.</p>
</div>
Built with <a href="https://www.sphinx-doc.org/">Sphinx</a> using a
<a href="https://github.com/readthedocs/sphinx_rtd_theme">theme</a>
provided by <a href="https://readthedocs.org">Read the Docs</a>.
</footer>
</div>
</div>
</section>
</div>
<script>
jQuery(function () {
SphinxRtdTheme.Navigation.enable(false);
});
</script>
</body>
</html>