these modifications is that modified EAM (MEAM) potentials pick up information
on the local order beyond the coordination number. For example, the nearest-
neighbor shell-tensor that can be associated with Eq. (62d) disappears for any
mirror-inversion crystal but not for the diamond lattice. Thus, MEAM-
potentials for silicon and germanium contain the square of the left-hand side
of Eq. (62d).

We now argue how a MEAM potential can stabilize open crystal lattices. To this
end we consider the expression

$\displaystyle g_{i}$ $\displaystyle=$ $\displaystyle{\sum_{j,k\neq
i,j}}\left\\{\cos\vartheta_{0}-n_{ij}^{\alpha}n_{ik}^{\alpha}\right\\}^{2}f_{j}(r_{ij})f_{k}(r_{ik})$
(63) $\displaystyle=$
$\displaystyle\rho_{i}^{2}\cos^{2}\vartheta_{0}-2\rho_{i\alpha}\rho_{i\alpha}\cos\vartheta_{0}+\rho_{i\alpha\beta}\rho_{i\alpha\beta}-(1-\cos\vartheta_{0})^{2}\sum_{j}f_{j}^{2}(r_{ij}),$

where the upper index on the partial embedding charge densities was omitted.
$g_{i}$ clearly has a minimum at the preferred bond angle $\vartheta_{0}$,
since $n_{ij}^{\alpha}n_{ik}^{\alpha}$ is the cosine of the angle $j$-$i$-$k$.
Now let us proceed with an ideal crystal structure, in which only partial
densities from nearest neighbors substantially contribute to the shell. In the
given approximation, we can replace $f(r_{ij})$ with $\rho/Z_{0}$ so that
after executing the restricted double sum correctly

$g_{i}\approx\rho^{2}\left(\cos^{2}\vartheta_{0}-2\cos\vartheta_{0}\nu_{0}^{\alpha}\nu_{0}^{\alpha}+\nu_{0}^{\alpha\beta}\nu_{0}^{\alpha\beta}\right)-\frac{\rho^{2}}{Z_{0}}(1-\cos\vartheta_{0})^{2}.$
(64)

Assigning an energy penalty on $g_{i}$ allows the degeneracy of different
crystalline structures with identical coordination number but different bond
angles to be lifted within a nearest-shell approximation. At the same time,
open structures with small $Z_{0}$ and correct bond angles can be favored over
closed-packed structures. It is clear that Eq. (64) can be constructed (with
the exception of the additive constant) from a linear combination of the MEAM
invariants like those given in Eq. (62).

Although MEAM potentials, similar in spirit to the ones just described, are
routinely used in molecular simulations, it is easily argued that these MEAMs
have issues. First, there is no physically motivated reason for the existence
of the upper indices in the $f_{j}^{(l)}$ coefficients when justifying the
modification from Baskes’ argument that the embedding function should also
depend on derivatives of the embedding density. Second, within this picture,
the embedding function can generally depend on scalars that can be constructed
from derivatives, which, however, must be constructed such that they are
invariant w.r.t. a rotation of the coordinate system. There are many more
scalars of a similar order as those contained in Eq. (62) [223]. For example,
if
$\rho_{\alpha_{1}...\alpha_{l}}(\mathbf{r})\coloneqq\partial_{\alpha_{1}}...\partial_{\alpha_{l}}\rho(\mathbf{r})$,
then in addition to $\rho_{\alpha\beta}^{2}$ there can be a contribution of
$\rho_{\alpha\alpha}\rho_{\beta\beta}$, which, in principle, can be linearly
independent of the former, rather than to be constrained to
$\rho_{\alpha\beta}^{2}-\rho_{\alpha\alpha}\rho_{\beta\beta}/3$, which can be
loosely associated with the term of Eq. (62c). Yet another invariant, which
can be constructed with four Cartesian indices and a square in the densities,
would be $\rho_{\alpha}\rho_{\alpha\beta\beta}$. Thus, to flesh out our second
point of criticism, the $l=2$ and $l=3$ terms in Eq. (62c) are far from
complete and it is not clear if the most relevant ones have been considered.

The “expansion” of Eq. (62c) truncates at the third-order term, which is the
lowest-order invariant that can discriminate between nearest-neighbor-shell
energies of fcc and hcp. However, different cubic structures cannot be
distinguished with shell tensors of rank three or less. Fourth, the
coefficients $t_{i}^{(l)}$ can depend, in principle, on $\rho_{i}^{(0)}$. It
might be possible to obtain some dependencies through asymptotic, dimensional,
or other systematic analyses. Guessing them correctly appears an almost
impossible task. Thus, we feel that the already impressive performance of
MEAMs should be further improvable. However, one potentially important
ingredient missing is the environment dependence of the functions
$f_{j}(r_{ij})$ beyond screening, e.g., an environment-dependence of
$\sigma_{\textrm{a}}$. Lifting this restriction would allow mediated
interaction to take place as they occur, for example, in breathing shell
models or in TB potentials going beyond the second-moment approximation,
albeit at a potentially much increased computational cost.

###  5.4 Stillinger-Weber potential and extensions

The binding energy in non-funky pair-, EAM, and TB2M potentials is largest
when atoms are most closely packed. However, the three lightest group-14
elements (C, Si, and Ge) are (meta-) stable in the diamond lattice, which has
a packing fraction less than 50% of that in fcc or hcp. To stabilize the
_open_ diamond lattice, Stillinger and Weber [228] proposed to add a rather
simple angle-dependent term to a Mie pair potential, which, however, was
multiplied with the rather clever cut-off function

$f_{\textrm{c}}(r_{ij},\gamma,r_{\textrm{c}})=\exp\left[\gamma/(r_{ij}-r_{\textrm{c}})\right]\,\Theta(r_{\textrm{c}}-r_{ij}).$
(65)

It has the nice property that all derivatives continuously approach zero at
the cut-off distance $r_{\textrm{c}}$. The angular add-on consisted of a
penalty quadratic in the deviation from tetrahedral bonding on a given atom
$i$, as originally proposed by Keating [229]. The extra-term is cut off when
the length of one of the two bonds, $r_{ij}$ and $r_{ik}$, forming a bond
angle of $\vartheta_{ijk}$ on atom $i$, approaches $r_{\textrm{c}}$. Put
together in the notation used throughout this article, the SW potential reads

$\begin{split}U(\\{\mathbf{r}\\})=&\sum_{i,j>i}\frac{U_{0}}{n-m}\left\\{m\left(\frac{r_{0}}{r_{ij}}\right)^{n}-n\left(\frac{r_{0}}{r_{ij}}\right)^{m}\right\\}f_{\textrm{c}}(r_{ij},\gamma_{1},r_{\textrm{c}})\\\
&+\sum_{i,j\neq i,k\neq
i,j}U_{\textrm{t}}\left\\{\cos\vartheta_{ijk}-\cos\vartheta_{\textrm{0}}\right\\}^{2}\,f_{\textrm{c}}(r_{ij},\gamma_{2},r_{\textrm{c}})\,f_{\textrm{c}}(r_{ik},\gamma_{2},r_{\textrm{c}})\end{split}$
(66)

with $\cos\vartheta_{\textrm{0}}=-1/3$ for the tetrahedral angle. Besides this
angle, the SW potential has eight independently adjustable parameters:
$U_{0}$, $r_{0}$, $m$, $n$, $\gamma_{1}$, $r_{\textrm{c}}$ for the pair
potential, and in addition $U_{\textrm{t}}$, $\gamma_{2}$ for the tetrahedral
bonding part. Although the SW potential is occasionally used to also model
carbon, it should be kept in mind that the original SW potential does not
recognize graphite as being energetically favorable at ambient pressure.

Despite its simplicity and empirical nature, the SW potential can be
parametrized for silicon to reproduce (by construction) its elastic properties
in the diamond structure, but also fairly well some high-pressure phases, and
the melting temperature at ambient pressure, including the anomalous jump from
small to large density during melting [228, 230, 231]. The latter goes hand in
hand with a coordination change from $Z_{0}=4$ in the diamond structure to
$Z\gtrsim 6$ in the liquid [232].

An appealing aspect of the SW potential is that the competition between
tetrahedral and dense packing favored by the angular and the pair part,
respectively, can be tuned through the ratio $U_{\textrm{t}}/U_{0}$ [233]. By
increasing it from a value characteristic for germanium through that of
silicon, it is possible to make the disordered phase adopt a local tetrahedral
order at large $U_{\textrm{t}}/U_{0}$, which explains why SW is occasionally
used for the united-atom modeling of water [234]. Using small
$U_{\textrm{t}}/U_{0}$, SW can describe condensed phases of tin, including
that condensing in the $\beta$-tin structure [235], which can be seen as a
compromise between tetrahedral and close packing.

Of course, as Stillinger and Weber [228] admitted, their model has
quantitative deficiencies. They were revealed most clearly by Biswas and
Hamann [236], who demonstrated that the equation of state of any crystalline
structure other than the diamond structure is highly flawed. However, even for
the diamond structure, the deficiency of the SW can become qualitative. For
example, SW fails to reproduce that silicon is brittle under tension, although
SW can be used to mimic brittle (non-silicon) tetrahedral solids, see Refs.
[237, 238]. This weakness is shared by other empirical potentials and related
to the finite range of the cutoff, which will be discussed in more detail in
Sect. 9.

The SW potential formally looks like a cluster potential. Given that the
three-body term has a zero-energy contribution in the ideal diamond structure,
it must be concluded that the pair interaction is an effective $Z_{0}=4$ pair
potential rather than the “true” pair potential. Improvements of the SW
address this issue partially by augmenting SW with an environment dependence
in terms of a coordination-dependent equilibrium bond angle $\vartheta_{0}$.
Such environment-dependent interaction potentials (EDIPs) were first developed
for silicon [239, 240] and later for other elements including carbon [241,
242]. Due to its ability to reflect the quasi-degeneracy of graphite and
diamond, an EDIP for carbon was the first empirical interatomic potential that
correctly predicted the relation of $3$-fold (graphite-like, sp2-hybridized)
and $4$-fold (diamond-like, sp3-hydridized) atoms in amorphous carbon obtained
from liquid quenches at varying densities [242].

Another interesting generalization of the SW potential, which Biswas and Haman
[236] proposed to model bond-angle energetics on atom $i$ beyond the Keating
model, reads

$U_{3}(r_{ij},r_{ik},r_{jk})=\sum_{l}C_{l}F_{l}(r_{ij},r_{ik})P_{l}(\cos\vartheta_{l}),$
(67)

where $C_{l}$ are expansion coefficients, $P_{l}(\cos\vartheta)$ Legendre
polynomials, and the $F_{l}(r_{ij},r_{ik})$ functions, which can be defined to
reproduce the exact or any general three-body potential. When assuming them to
factorize as $F_{l}(r_{ij},r_{ik})=\phi_{l}(r_{ij})\phi_{l}(r_{ik})$, e.g.,
with simple exponential dependencies, Biswas and Hamann demonstrated the
computation of $U_{3}$ to be reducible to order $Z_{\textrm{loc}}$, where
$Z_{\textrm{loc}}$ is the number of atoms within the interaction range of atom
$i$, while irreducible potentials, such as the ATM potential, require a number
of operations proportional to $Z_{\textrm{loc}}^{2}$.

###  5.5 Tersoff potentials

Dissatisfied by the poor transferability of the SW potential [236] and
inspired by Abell’s analysis of the sensitivity of bond strengths on the local
environment [62], Tersoff presented a series of papers [243, 244, 245, 246,
247], in which he introduced the concept of bond order to the world of
empirical potentials. According to Pauling [248], the bond order is the
difference of the number of electrons in bonding and anti-bonding orbitals. It
is a monotonically decreasing function of the coordination number. Tersoff’s
work was an attempt to construct functions that measure (effective)
coordination numbers and thereby the bond order. To do so, he introduced a
cut-off function, which is unity up to a distance $r_{\textrm{c}1}$ meant to
include nearest neighbors and which then quickly falls off to zero at a
distance $r_{\textrm{c}2}$ supposedly less than typical next-nearest neighbor
distances:

$f_{\textrm{c}}(r)=\Theta(r_{\textrm{c}1}-r)+\frac{\Theta(r-r_{\textrm{c}1})\Theta(r_{\textrm{c}2}-r)}{2}\left\\{1+\cos\left(\pi\frac{r-r_{\textrm{c}1}}{r_{\textrm{c}2}-r_{\textrm{c}1}}\right)\right\\}.$
(68)

With the help of this cut-off function, the bonds to atom $i$ other than the
$ij$ bond can be characterized with

$\xi_{ij}=\sum_{k\neq i,j}g(\cos\vartheta_{ijk})f_{\textrm{c}}(r_{ik}),$ (69)

where the angular dependence, originally chosen as

$g(\cos\vartheta)=1+\frac{c^{2}}{d^{2}}-\frac{c^{2}}{d^{2}+(\cos\vartheta_{0}-\cos\vartheta)^{2}}$
(70)

will later help to make a crystal in the diamond structure resist shear
stresses. Here $c$ and $d$ are adjustable parameters, while $\vartheta_{0}$
is, for example, the ideal tetrahedral angle. The term $\xi_{ij}$ can be
interpreted as an effective coordination number (minus one, since atom $j$ is
excluded from the sum) under which atom $i$ “sees” atom $j$. If all bond
angles on atom $i$ are equal to $\vartheta_{0}$, $\xi_{ij}=Z_{i}-1$ assuming
the cut-off function to assume the values one and zero for nearest and next-
nearest neighbors, respectively. However, “unhappy” angles formed by $ij$ and
$ik$ bonds lead to an increase of $\xi_{ij}$.

Tersoff then constructed a prefactor, $b_{ij}$, to the attractive pair
interaction as a function of the (effective) coordination number and called it
a bond-order variable. It is defined heuristically as

$b_{ij}(\xi_{ij})=\left\\{1+(\beta_{i}\xi_{ij})^{\gamma_{i}}\right\\}^{-1/(2\gamma_{i})},$
(71)

where $\beta_{i}$, and $\gamma_{i}$ are element-specific parameters. For large
$\xi_{ij}$, $b_{ij}$ scales with $1/\sqrt{\xi_{ij}}$, which was shown by
Brenner [249] to be similar to EAM or TB2M potentials. With these ingredients,
a general Tersoff potential reads

$U(\\{\mathbf{r}\\})=\frac{1}{2}\sum_{i,j\neq
i}\left\\{u_{\textrm{r}}(r_{ij})-b_{ij}u_{\textrm{a}}(r_{ij})\right\\}f_{\textrm{c}}(r_{ij}),$
(72)

where the (pair) contributions $u_{\textrm{r}}(r_{ij})$ and
$u_{\textrm{a}}(r_{ij})$ are strictly repulsive and strictly attractive,
respectively. Typically, $u_{\textrm{a,r}}(r_{ij})$ are exponential functions
to reflect the universal equation of state of crystals [250, 62].

Using a symmetrized bond-order parameter, $\bar{b}_{ij}=(b_{ij}+b_{ji})/2$,
the Tersoff potential can be cast as

$U(\\{\mathbf{r}\\})=\sum_{i,j>i}\frac{U_{0}}{n-m}\left\\{m\,e^{n(1-r_{ij}/r_{0})}-\bar{b}_{ij}\,n\,e^{m(1-r_{ij}/r_{0})}\right\\}f_{\textrm{c}}(r_{ij}).$
(73)

Tersoff chose $n=2m$, as to reproduce the original Morse potential [246]. From
the Ducastelle-potential perspective, i.e., when setting $g(\cos\vartheta)$ to
unity, this is precisely the choice of exponents, in which the binding energy
is insensitive to the coordination number within the nearest-neighbor
approximation, see our discussion on this in Sect. 6.1. Tweaking the exponents
toward larger (smaller) $n/m$ ratios biases structures in favor of larger
(smaller) coordination numbers.

A point that might be particularly important and easy to make in the context
of bond-order potentials is that the effect of screening is cast through a
cutoff rather than through the analysis of local topology. For example, given
fixed values for the two cut-off radii $r_{\textrm{c}1,2}$ in the Tersoff
potential, it simply seems wrong that the variable $\xi_{ij}$ can change from
its generic value of 3 in the diamond structure to 15 upon a relatively minor
isotropic compression or to 0 upon a hypothetical homogeneous decompression. A
more elaborate discussion of this issue is given in Sect. 9.

###  5.6 Generic functional form

All open-shell potentials discussed in the preceding sections, and the closed-
shell ATM potential, can be cast into the universal functional form given by
Eqs. (53) and (54). The universal form highlights similarities between the
construction of these empirical potentials and simplifies the analytic
manipulation, such as the computation of properties that depend on derivatives
like forces,

$\displaystyle\mathbf{f}_{n}=$ $\displaystyle-$
$\displaystyle\sum_{j}\left(\frac{\partial U_{2}}{\partial
r_{nj}}+\frac{\partial U_{\text{m}}}{\partial
r_{nj}}\right)\frac{\mathbf{r}_{ij}}{r_{ij}}$ (74) $\displaystyle-$
$\displaystyle\frac{1}{2}\sum_{ij}\frac{\partial
U_{\text{m}}}{\partial\xi_{ij}}\sum_{k\not=i,j}\left\\{\mathbf{\Xi}_{ij}(\delta_{in}-\delta_{jn})+\mathbf{\Xi}_{ik}(\delta_{in}-\delta_{kn})+\mathbf{\Xi}_{jk}(\delta_{jn}-\delta_{kn})\right\\}$

with $\mathbf{\Xi}_{ij}\equiv(\partial\Xi/\partial r_{ij})\,\hat{n}_{ij}$,
where $r_{ij},r_{ik},r_{jk}$ are the first, second, and third argument of
$\Xi(...)$. The remaining algebraic operations are partial derivatives in
scalar variables that are straightforward to carry out. Stresses or higher
order derivatives like elastic constants are equally straightforward to
evaluate.

We summarize the functional forms for the potentials discussed here in Tab. 3.
Both ATM and SW potentials have a many-body contribution of the form
$U_{\text{m}}(r_{ij},\xi_{ij})=\xi_{ij}$, which reduces the many-body
contribution to a true three-body contribution, e.g., $U_{3}$ in Eq. (2).
Carlsson [19] dubbed these types of expressions _cluster potentials_.
Ducastelle and Tersoff introduce a nonlinear mapping for the functional
dependency on $\xi_{ij}$ in $U_{\text{m}}$, turning the cluster potential into
what Carlsson termed a _cluster functional_. The critical difference is that
while cluster potentials go to a finite order in the slowly converging series
Eq. (2), cluster functionals implicitly include higher-order terms. This can
be seen by expressing the many-body contribution as a Taylor series in
$\xi_{ij}$. This Taylor series is truncated for cluster potentials but
contains quadratic $\xi_{ij}^{2}$, cubic $\xi_{ij}^{3}$ and higher-order terms
for cluster functionals. The quadratic term,

$\xi_{ij}^{2}=\sum_{k,l}\Xi(r_{ij},r_{ik},r_{jk})\Xi(r_{ij},r_{il},r_{jl}),$
(75)

clearly contains a four-body contribution to $U_{4}$ in Eq. (2). Similarly,
the cubic terms contributes a five-body interaction and so on. This implicit
incorporation of higher-order terms can alleviate the slow convergence of the
formal series expansion given by Eq. (2), while avoiding computationally
expensive, explicit calculations of higher-order terms.

$\displaystyle U(\\{\mathbf{r}\\})$
$\displaystyle=\frac{1}{2}\sum_{ij}\left\\{U_{2}(r_{ij})+U_{\text{m}}(r_{ij},\xi_{ij})\right\\}\quad\text{with}\quad\xi_{ij}=\sum_{k\not=i,j}\Xi(r_{ij},r_{ik},r_{jk})$
$\displaystyle U_{2}(r_{ij})$ $\displaystyle\propto\begin{cases}\text{$m$-$n$
Lennard-Jones}&\text{ATM, SW}\\\
\exp(-r_{ij}/\sigma_{\text{r}})f_{\text{c}}(r_{ij})&\text{Ducastelle/Gupta,
Tersoff}\end{cases}$ $\displaystyle U_{\text{m}}(r_{ij},\xi_{ij})$
$\displaystyle\propto\begin{cases}\xi_{ij}&\text{ATM, SW}\\\
\delta_{ij}\sqrt{\xi_{ij}}&\text{Ducastelle/Gupta}\\\
(1+\xi_{ij})^{-1/(2\gamma)}\exp(-r_{ij}/\sigma_{\text{a}})f_{\textrm{c}}(r_{ij})&\text{Tersoff}\end{cases}$
$\displaystyle\Xi(r_{ij},r_{ik},r_{jk})$
$\displaystyle\propto\begin{cases}(1+3\cos\vartheta_{ijk}\cos\vartheta_{jki}\cos\vartheta_{kij})/(r^{3}_{ij}r^{3}_{jk}r^{3}_{ki})&\text{ATM}\\\
\exp(-r_{ik}/\sigma_{\text{a}})&\text{Ducastelle/Gupta}\\\
g^{\text{SW}}(\cos\vartheta_{ijk})f^{\text{SW}}_{\textrm{c}}(r_{ij})f^{\text{SW}}_{\textrm{c}}(r_{ik})&\text{SW}\\\
g^{\text{T}}(\cos\vartheta_{ijk})f_{\textrm{c}}(r_{ik})&\text{Tersoff}\end{cases}$
Table 3:  Selected many-body potentials as cast into the generic functional
form, Eqs. (53) and (54). We here present simplified expressions for unary
systems and have omitted prefactors to highlight functional dependencies, but
generalization to alloys and introduction of prefactors is straightforward.
$\sigma_{\text{r}}$ and $\sigma_{\text{a}}$ are characteristic length-scales
of repulsion and attraction, respectively. $f_{\text{c}}$ are cutoff functions
whose expressions differ between Stillinger-Weber and Tersoff, see Eq. (131).
The functions $g^{\text{SW}}$ and $g^{\text{T}}$ encode the angular dependence
of Stillinger-Weber and Tersoff, respectively. Divergent on-site terms
$\xi_{ii}$ are implicitly excluded from summation. Generic forms of the SW and
Tersoff potentials differ from their original formulation by a (constant)
onsite contribution that is inconsequential for forces.

###  5.7 Brenner potentials

Brenner reparameterized Tersoff’s functional form, specifically with
$\gamma_{i}=1$, for modeling hydrocarbon chemistry [10, 251]. This potential
is now known as the _reactive empirical bond-order potential_ (REBO). Brenner
realized that Tersoff’s potential essentially ignores the $\pi$ bond in
$sp$-valent system. Thus, in its raw form, it is unable to describe the
complex chemistry of hydrocarbons, for example radicals or bond conjugation.
To correct for these deficiencies, Brenner introduced lookup tables that
correct Tersoff’s total energy expression depending on the coordination of the
atoms involved in the bond. The many-body term becomes

$U_{\text{m}}(r_{ij},\xi_{ij})=-\left[b_{ij}(\xi_{ij})+\Pi_{ij}(N_{i},N_{j},N_{ij}^{\text{conj}})\right]u_{\text{a}}(r_{ij})f_{\text{c}}(r_{ij}),$
(76)

where $N_{i}$ is the coordination of atom $i$ obtained with a smooth cutoff
function $f_{\text{c}}(r)$, $N_{i}=\sum_{j}f_{\text{c}}(r_{ij})$. The function
$\Pi_{ij}$ accounts for the missing $\pi$-bond and is a bicubic interpolation
between values tabulated for specific coordination of the atoms.
$N_{ij}^{\text{conj}}$ is another integer variable that is $\geq 2$ if the
bond is conjugated, which is the case if a neighboring atom (of the bond
$i$-$j$) has a coordination less than four. $N_{ij}^{\text{conj}}$ increases
the neighbor shell, which influences the energy of bond $i$-$j$, making the
potential longer ranged than Tersoff’s.

In addition to the correction for the $\pi$-bond, Brenner also corrected the
$\sigma$-bond contribution, using the bond-order expression

$b_{ij}(\zeta_{ij})=\left[1+\xi_{ij}+P_{ij}(N_{i}^{\mathrm{(H)}},N_{i}^{\mathrm{(C)}})\right]^{-1/2},$
(77)

where $P_{ij}$ is again a lookup table and $N_{i}^{\mathrm{(H)}}$ and
$N_{i}^{\mathrm{(C)}}$ are coordination numbers of atom $i$, but only counting
hydrogens or carbons, respectively, such that
$N_{i}=N_{i}^{\mathrm{(H)}}+N_{i}^{\mathrm{(C)}}$. $P_{ij}$ is only nonzero
for C-C bonds if a hydrogen is present as a bonding partner,
$N_{i}^{\mathrm{(H)}}>0$. The lookup tables were fitted to atomization
energies of many different hydrocarbon molecules. More information can be
found in Refs. [10, 251], including nodal values for the lookup table
$\Pi_{ij}$ and $P_{ij}$.

A _second generation reactive empirical bond-order potential_ (REBO2) was
published twelve years later by Brenner et al. [11]. This potential improves
upon the original formulation by adopting new functional forms for
$u_{\text{r}}(r)$ and $u_{\text{a}}(r)$. Specifically, the repulsion in REBO2
contains a screened Coulomb (Yukawa) contribution, as discussed in Sect.
3.2.3, in addition to the purely exponential repulsion of REBO and Tersoff’s
formulation. REBO2 has two distinct angular functions $g(\vartheta_{ijk})$,
depending on the coordination of atom $i$. The potential also added a (four-
atom) dihedral potential, allowing proper modeling of rotations around carbon-
carbon double bonds. The REBO2 potential was augmented with non-bonded
(dispersion) interactions by Stuart, Tutein and Harrison [252] (a development
that, despite being published earlier, occurred after the development of
REBO2).

We here put forward the bold claim that REBO is the first successful machine-
learned potential. Essentially, Brenner did away with the fixed functional
form through his lookup tables to correct Tersoff’s potential for a specific
application. While the sophisticated interpolation techniques that are
nowadays used for machine-learned potentials (neural networks, Gaussian
processes [253]) were still under development at the time of Brenner’s
publications, the spline-interpolated lookup tables essentially fulfill a
similar role of extrapolating in a high-dimensional space.

###  5.8 Beyond second moments

Before talking about potentials _beyond_ second moments, we need to take a
step back and elucidate the origins of the second-moment expansion. In
particular, we want to emphasize that there is an atom-centered and a bond-
centered expansion that lead to different functional forms for binding
energies. We start our discussion with the atom-centered expansion.

The TB2M expressions can be rationalized from orthogonal tight-binding models
as follows: The total (band) energy of the electronic system is given by

$E_{\text{band}}=\tr\hat{\rho}\hat{H}=\sum_{n}f(\varepsilon_{n})\varepsilon_{n}=\rho_{i\alpha
j\beta}H_{i\alpha j\beta}$ (78)

where the $f(\varepsilon)$ is the Fermi-Dirac distribution, $\varepsilon_{n}$
are the energy eigenvalues and $\hat{\rho}=f(\hat{H})$ is the density
matrix.222Functions of the Hamitonian are defined through their action on
eigenfunctions, $f(\hat{H})|n\rangle=f(\varepsilon_{n})|n\rangle$, or through
a power series expansion of $f$. At low electronic temperature,
$f(\varepsilon)$ is simply a step function such that the sum in Eq. (78) runs
over all states with energy below the Fermi level $\varepsilon_{\text{F}}$.
The total energy of the system typically includes a pair-wise repulsive
contribution in addition to $E_{\text{band}}$, which absorbs (and
approximates) all the things that were implicity or explicitly neglected in
the derivation of the orthogonal tight-binding band model, such as three-
center terms, overlap repulsion, and others.

A crucial step in turning the tight-binding total energy picture into the atom
or bond-centric pictures underlying interatomic potential was the formulation
of the local density of states (DOS) based on moments of the Hamiltonian. In
metals, bonding is determined by the shape of the electronic density of states
around the Fermi level. The local density of states is given by [34]

$n_{i\alpha}(\varepsilon)=\langle
i\alpha|\delta(\varepsilon-\hat{H})|i\alpha\rangle,$ (79)

where $|i\alpha\rangle$ is orbital $\alpha$ on atom $i$. In terms of this
local density of states, the band energy becomes

$E_{\text{band}}=\sum_{i\alpha}\int\operatorname{d\\!}\varepsilon\,\varepsilon
f(\varepsilon)n_{i\alpha}(\varepsilon).$ (80)

The _moments_ of the local density of states are given by [254, 255],

$\mu_{i\alpha}^{(p)}=\int\operatorname{d\\!}\varepsilon\,\varepsilon^{p}n_{i\alpha}(\varepsilon)=\langle
i\alpha|\hat{H}^{p}|i\alpha\rangle=H_{i\alpha j\beta}H_{j\beta k\gamma}\cdots
H_{l\delta i\alpha},$ (81)

without a sum over $i\alpha$ (although that index is repeated). The right hand
side of Eq. (81) has an intuitive interpretation. For the second moments,

$\mu_{i\alpha}^{(2)}=H_{i\alpha j\beta}H_{j\beta i\alpha},\qquad\text{(no
summation over $i\alpha$!)}$ (82)

we hop from atom $i$ to all neighboring atoms $j$ and back to atom $i$ and the
sum of the squared Hamiltonian elements given the moment. For the third
moment,

$\mu_{i\alpha}^{(3)}=H_{i\alpha j\beta}H_{j\beta k\gamma}H_{k\gamma i\alpha},$
(83)

we take all self-returning hopping paths involving three atoms, and so on.
Since the Hamiltonian matrix is short ranged (remember, it decays roughly like
an exponential!), the second moment only depends on the close vicinity of atom
$i$. Higher moments “feel” the atomic structure out to further distances. The
moments are useful because they characterize the shape of the density of
states. If we know all moments, we know the density of states; conversely we
can approximate the density of states using just the lowest moments.

Systematically approximating the density of states from the lowest moments is
highly nontrivial, and a number of approaches have been suggested. These
approaches include the recursion method [256, 257] (that lead to the
systematic development of _bond-order_ potentials [258, 259]), the maximum
entropy method [260, 261, 262, 263], or simply approximating the local DOS
with a rectangle [264, 204]. Let us now assume we have an unary (for
simplicity $s$-valent) system (e.g. Cu) with bond integral
$\beta_{ss\sigma}(r)$ for the ss$\sigma$-bond. There is no angular dependence
and the second moment of the local DOS at atom $i$ – see Eq. (82) – is given
by

$\mu^{(2)}_{is}=\sum_{j}\beta_{ss\sigma}^{2}(r_{ij}).$ (84)

Assuming a half-filled rectangular DOS and a single electron per atom in our
$s$-valent metal, the contribution of atom $i$ to the band _energy_ (the first
moment of the local DOS integrated up to the Fermi level) becomes

$E_{\text{band},i}=-\sqrt{6\mu_{is}^{(2)}},$ (85)

which is the origin of the square-root in Eq. (59a). A maximum entropy
estimate would lead to a Gaussian, which changes the numerical prefactor of
Eq. (85) but not the functional form [265].

Semiempirical atom-centered higher-moment methods have also been developed.
Carlsson considered fourth moments (that define something like the kurtosis of
the DOS) for models of semiconductors [266] and transition metals [267]. While
higher-moment atom-centered techniques have seen some development since [268,
269], they have to the best of our knowledge not been widely employed in
molecular calculations. The breakthrough of higher-moment methods came with
the development of systematic bond-centered techniques, briefly described in
the following.

The arguments present above imply an atom-centered view, as the local density
of states is a per-atom quantity. For a slight shift in perspective, we
abandon this band-centric view and define
$E_{\text{band}}=E_{\text{bond}}+E_{\text{prom}}$ where the _bond energy_

$E_{\text{bond}}=\frac{1}{2}\sum_{i\not=j}2\rho_{i\alpha j\beta}H_{i\alpha
j\beta}$ (86)

and the _promotion energy_ $E_{\text{prom}}$ contains the diagonal terms, i.e.
the on-site energies of the tight-binding model. This construction is the
basis of the tight-binding _bond_ model [270]. Equation (86) for the bond
energy is appealing, because it introduces the bond $i$-$j$ as the central
quantity: It has the appearance of a pair interaction with pair energy
$E_{ij}=2\rho_{i\alpha j\beta}H_{i\alpha j\beta}$ – compare it for example to
the attractive part of the Tersoff potential given by Eq. (72). The term
$2\rho_{i\alpha j\beta}$ is called the _bond order_. As we know from Tersoff’s
potential, this bond order – or rather the bond order variable $b_{ij}$, Eq.
(71) – depends on the environment of the specific bond and modulates its
strength. The promotion energy is constant if local charge neutrality is
imposed and can be ignored for the computation of forces and stresses [259].

We can write the density matrix as

$\rho_{i\alpha
j\beta}=\int\operatorname{d\\!}\varepsilon\,f(\varepsilon)n_{i\alpha
j\beta}(\varepsilon)\quad\text{with}\quad n_{i\alpha
j\beta}(\varepsilon)=\langle
i\alpha|\delta(\varepsilon-\hat{H})|j\beta\rangle,$ (87)

which is easily seen by inserting this expression into Eq. (78). Note that
written this way, the density matrix has a form similar to Eqs. (79) and (80).
The onsite terms $\rho_{i\alpha i\alpha}$ correspond to the integral over the
local density of states $n_{i\alpha}(\varepsilon)\equiv n_{i\alpha
i\alpha}(\varepsilon)$ up to the Fermi level, which are the number of
electrons $N_{i\alpha}$ in that orbital. We now consider the linear
combination of orbitals
$|\pm\rangle=(|i\alpha\rangle\pm|j\beta\rangle)/\sqrt{2}$ that correspond to
bonding ($+$) and anti-bonding ($-$) states [205]. Then
$n_{\pm}=\frac{1}{2}\left\\{n_{i\alpha}+n_{j\beta}\right\\}\pm n_{i\alpha
j\beta}$ and hence $n_{i\alpha j\beta}=(n_{+}-n_{-})/2$. The bond-order can
therefore be written as

$2\rho_{i\alpha j\beta}=N_{+}-N_{-},$ (88)

with $N_{\pm}=\int\operatorname{d\\!}\varepsilon
f(\varepsilon)n_{\pm}(\varepsilon)$. This is the difference of electrons in
anti-bonding ($N_{-}$) and bonding ($N_{+}$) orbitals. We conclude that
Pauling’s interpretation [248] of the bond-order was indeed correct.

Exploiting the symmetry of the $i$-$j$ bond, the off-diagonal Hamiltonian
block $H_{ij}$ can be diagonalized (by rotating into the frame of the local
bond $i$-$j$), allowing the pair energy to be written as [271],

$E_{ij}=2H_{\sigma}(r_{ij})\rho_{\sigma,ij}+2H_{\pi}(r_{ij})\rho_{\pi,ij}+2H_{\delta}(r_{ij})\rho_{\delta,ij},$
(89)

the sum of $\sigma$, $\pi$ and $\delta$ contributions to the bond energy.
Pettifor, Aoki, Horsfield and coworkers [258, 272, 273, 259] developed a
systematic expansion of the bond-order in the moments computed on atom $i$ and
$j$. To second order, they find
$\rho_{\sigma,ij}\propto(\mu_{i\sigma}^{(2)}+\mu_{j\sigma}^{(2)})^{-1/2}$,
leading to an expression similar to Tersoff’s total energy (if $\pi$-bonds are
neglected). Since this pioneering work, higher-order bond-order potentials
have been developed for a number of elements and elemental combinations, such
as hydrocarbons [274, 275, 276], molybdenum [277], gallium-arsenide [278],
silicon [279], tungsten [280], iron [281], cadmiun-telleride [282], and
aluminum-copper-hydrogen [283, 284, 285].

##  6 Properties of pair and many-body potentials

###  6.1 Binding energy, lattice constants, and equation of state

The properties of crystals with a single atom in the primitive cell that
interact through a smooth, short-range two-body potential $U_{2}(r)$ [with a
single (steep) minimum in $U_{2}(r)$ and a single extremum in relevant higher-
order derivatives] are severely restricted. We start by discussing the
cohesive energy, arguably the most important property of any structure. For a
perfect mono-atomic and thus non-ionic crystal, in which each atom is
equivalent, the energy per atom in the (effective) pair-potential
approximation can be given by

$U_{\textrm{pa}}(a_{0})=\frac{1}{2}\sum_{s=0}Z_{s}U_{2}(a_{s}).$ (90)

Values for $Z_{0}$, $Z_{1}$, $a_{1}/a_{0}$ are stated in Table 2 for mono-
atomic cubic crystals along with other characteristic, dimensionless numbers,
which will be introduced further below. Given that $Z_{0}=12$ is the maximal
coordination number, face-centered cubic (fcc) and hexagonal closed packed
(hcp), which has the same $Z_{0}$, $Z_{1}$, and $a_{1}/a_{0}$ as fcc, are the
generally preferred structures of mono-atomic systems with short-range
interactions. What phase wins depends on the precise functional form of the
potential but also where and how the potentials are cut off, as discussed in
more detail in Sect. 9.1.

To lowest order, i.e., when only including interactions within the first
shell, the nearest-neighbor spacing satisfies $a_{0}=r_{0}$, where $r_{0}$ is
the location of the energy minimum of the dimer. Like nearest neighbors, more
distant neighbor shells are located at a negative potential energy. As a
result, the total binding energy per atom in the solid will exceed $Z_{0}$
times the energy per atom in the dimer, which is $U_{0}/2$. Due to the further
shells sitting in the attractive part of the potential, the lattice will
contract. As a consequence of this contraction, nearest neighbors are pushed
more deeply into their repulsive interaction, which then leads to an increase
of the bulk modulus $B$, despite counteracting effects of other shells. The
inequalities

$\displaystyle a_{0}^{\textrm{eq}}/r_{0}$ $\displaystyle\lesssim 1$ (91a)
$\displaystyle U_{\textrm{pa}}^{\textrm{eq}}(Z_{0})/U_{0}$
$\displaystyle\gtrsim Z_{0}/2$ (91b) $\displaystyle B^{\textrm{eq}}$
$\displaystyle\gtrsim{Z_{0}}\,\frac{k_{\textrm{d}}r_{0}^{2}}{18\,v_{\textrm{pa}}^{\textrm{eq}}}$
(91c)

encapsulate this discussion, $U_{\textrm{pa}}^{\textrm{eq}}$ being the binding
energy per atom (in the ground state with nearest-neighbor distance
$a_{0}^{\textrm{eq}}$) and $v_{\textrm{pa}}^{\textrm{eq}}$ is the volume per
atom. Values of $v_{\textrm{pa}}$ for different structures are tabulated in
Tab. 2. These equations are expected to hold for non-ionic solids interacting
through pair potentials. $k_{\text{d}}$ is the curvature of the pair potential
in the minimum, which could be called the dimer spring constant. It takes the
value

$k_{\textrm{d}}=mnU_{0}/r_{0}^{2}$ (92)

for the pair potentials considered here, at least in their most simple form
without amendments. The presented treatment is easily generalized to simple
binaries as long as Coulomb interactions are negligible.

Ionic crystals with long-range Coulomb interaction must be treated separately,
because $Z_{s}$ can be of order $Z_{0}a_{s}/a_{0}$ so that the sum over shells
is not convergent. Summing up Coulomb interactions (in periodically repeated
systems) is an entire industry by itself [286, 287]. However, as long as
crystallographic positions are fixed relative to the unit cell, it is clear
that the sum over $1/r^{m}$ terms is a geometric factor depending on the
lattice times $1/a_{0}^{m}$. To keep the treatment analytically amenable, we
only consider nearest-neighbor repulsion plus Coulomb interaction and
approximate the energy per atom as

$U_{\textrm{pa}}(a_{0})\approx\frac{Z_{0}}{2}\,\frac{U_{0}}{n-1}\,e^{n(1-a_{0}/r_{0})}-\frac{\alpha_{\textrm{M}}}{8\pi\varepsilon_{0}}\,\frac{Q^{2}}{a_{0}},$
(93)

where $\alpha_{\textrm{M}}$ is the geometric factor for $m=1$, the Madelung
constant. Since $\alpha_{\textrm{M}}<Z_{0}$ for all (simple) crystal
structures, it is clear that the absolute ratio of Coulomb and repulsive
potential is reduced compared to that of the dimer. This leads to a reduction
of the energy per atom in the crystal and consequently to an inversion of all
three inequalities listed in Eq. (91a). Eq. (93) is not necessarily
sufficiently accurate to determine what phase an ionic crystal assumes as a
function of the parameter $n$ [152, 48]. To this end, it is necessary to also
include dispersive interaction and repulsion between anions plus potentially
many-body terms [152], which are frequently cast through so-called breathing
and deformable shell models as described in Sect. 4.4.

Using Mie instead of exponential repulsion and restricting repulsion to
nearest neighbors allows simple closed-form expressions to be derived [155],

$\displaystyle a_{0}^{\textrm{eq}}/r_{0}=$
$\displaystyle\sqrt[n-1]{Z_{0}/\alpha_{\textrm{M}}}$ (94a) $\displaystyle
U_{\textrm{pa}}^{\textrm{eq}}/U_{0}=$
$\displaystyle\sqrt[n-1]{\alpha_{\textrm{M}}/Z_{0}}\,\alpha_{\textrm{M}}$
(94b) $\displaystyle B^{\textrm{eq}}=$
$\displaystyle\sqrt[n-1]{\alpha_{\textrm{M}}/Z_{0}}\,\alpha_{\textrm{M}}\,\frac{n}{n-1}\,\frac{k_{\textrm{d}}r_{0}^{2}}{18\,v_{\textrm{pa}}},$
(94c)

where $k_{\textrm{d}}$ follows from Eq. (92). As expected, the gain in energy
and stiffness obtained through the condensation from molecules into crystals
is much reduced compared to that of short-range potentials.

Given $U_{\textrm{pa}}(a_{0})$, the equation of state (EOS) at small
temperature follows from $p=-\partial U_{\textrm{pa}}/\partial
v_{\textrm{pa}}$, which can be evaluated quite easily through

$p(v_{\textrm{pa}})=-\frac{\partial U_{\textrm{pa}}}{\partial
a_{0}}\,\frac{\partial a_{0}}{\partial v_{\textrm{pa}}}$ (95)

when using the volumes per atom $v_{\textrm{pa}}(a_{0})$ listed in Tab. 2. The
lean pair potentials introduced here allow the EOS of our four reference
compounds to be reproduced quite accurately in their thermodynamically stable
phase, as is revealed in Fig. 10. However, for copper and carbon, the obtained
parameters differ quite substantially from those describing the dimers, while
those for NaCl and Ar differ only by roughly 10–20%. Pertinent numbers are
listed in Tab. 1.

Figure 10: Equation of state, $p(V)$, for C, Cu, Ar, and NaCl as determined
experimentally and by pertinent fits assuming pair-additivity of potentials of
the same functional form as in Fig. 3 and with parameters listed in Table 1.
The pressure $p$ and the volume $V$ are normalized to the zero-pressure bulk
modulus $B$ and volume $V_{0}$, respectively, which the fits were constrained
to reproduce. Black full circes [288] and open circles [289] show reference
data for fcc argon, orange plus signs [290] and crosses [291] for fcc copper,
red diamonds [292] for carbon in the cubic-diamond structure as well as blue,
full [293] and open [294] squares for NaCl in the rock salt structure.

The apparent dimer binding energies $U_{0}$ deduced from the fits are
decreased by a factor of approximately 2.2 for carbon ($Z_{0}=4$) and by a
factor of 4.6 for copper ($Z_{0}=12$). Both numbers can be crudely
rationalized using the EAM/TB2M potential proposed by Ducastelle (see Sect.
5.2), which predicts the binding energy per atom to increase a little less
slowly than with $\sqrt{Z_{0}}$.

Frequently, a measured or computed EOS is also fitted to the so-called
“universal equation of state” [250, 295, 296] given by

$p(V)=3\,B_{0}\,\frac{1-\eta}{\eta^{2}}e^{3\,(B_{0}^{\prime}-1)\,(1-\eta)/2}.$
(96)

Here, $B_{0}$ is the bulk modulus at $p=0$ and $\eta=\sqrt[3]{V/V(p=0)}$,
while $B^{\prime}_{0}$ is the change of the bulk modulus $\partial B/\partial
p$ at $p=0$ and constant temperature, which is ideally zero temperature when
gauging interaction parameters from lattice sums. In fact, the reference data
for carbon [292] in Fig. 10 was produced using Eq. (96) whose adjustable
coefficients had been gauged on experimental data.

As can be seen, also an ionic solid satisfies the universal EOS under
compressive stress reasonably well, although its two-body potential is
markedly distinct from that of carbon or copper. Thus, inverting an EOS, in
particular when obtained solely under compression, allows potentials to be
tested but not to be parametrized. Using first-principle methods, the EOS can
be extended to tensile stresses, thereby providing useful information also on
the attractive part of the potential. Careful parameterizations of potentials
make use of such data [297, 298, 299, 300].

Even though pair potentials can be used to fit the EOS of individual
crystalline structures, they are not transferable between them. As described
in detail in Sect. 5, many-body terms are required for transferability. We
will now analyze the Ducastelle potential as one of the simplest many-body
formulations for metals. In the nearest-neighbor approximation, the per-bond
energy of crystalline Ducastellium, Eq. (60), can be approximated as

$U_{\textrm{pb}}(a_{0},Z_{0})=\frac{U_{0}}{m-n}\left\\{me^{n(1-a_{0}/r_{0})}-\frac{1}{\sqrt{Z_{0}}}ne^{m(1-a_{0}/r_{0})}\right\\}.$
(97)

It follows for the equilibrium bond length that

$\frac{a_{0}^{\textrm{eq}}(Z_{0})}{r_{0}}=1+\frac{1}{2}\frac{\ln Z_{0}}{n-m}$
(98)

and for binding energy

$U_{\textrm{pa}}^{\text{eq}}(Z_{0})=U_{0}Z_{0}^{\mu}/2\quad\text{or}\quad
U_{\textrm{pb}}^{\text{eq}}(Z_{0})=U_{0}Z_{0}^{\mu-1}$ (99)

with $\mu=(n/2-m)/(n-m)<1/2$. Moreover, the generalized stiffness of the
potential turns out to be

$k^{\text{eq}}(Z_{0})=Z_{0}^{\mu-1}k^{\text{eq}}(1),$ (100)

where $k^{\text{eq}}(1)$ is the stiffness of the bond in the dimer.
Corrections due to next-nearest and more distant shells alter these results.
However, all trends are robust and reflect the results for copper shown in
Fig. 9—as well as for a few other metals condensing in fcc. These trends are a
logarithmic increase of $a_{0}^{\text{eq}}$ with $Z_{0}$ and an algebraic
increase of $U_{\text{pa}}^{\text{eq}}$ with a power less than $1/2$. It is
important to note that for $n=2m$ we obtain $\mu=0$, and Eq. (99) shows that
the cohesive energy is then independent of crystal structure.

Note that the same trends are also observed by carbon (also shown in Fig. 9),
despite the fact that it cannot be described by Ducastelle’s potential, as
carbon’s ground-state structures are open crystal lattices. To understand
covalent bonding, it is instructive to consider the energy per bond

$U_{\text{pb}}(a_{0})=\frac{U_{0}}{n-m}\left\\{m\,e^{n(1-a_{0}/r_{0})}-\bar{b}\,n\,e^{m(1-a_{0}/r_{0})}\right\\}f_{\textrm{c}}(r_{ij})$
(101)

of nearest-neighbor Tersoffium, Eq. (73). We directly see that we recover Eq.
(97) of Ducastellium for $\bar{b}=1/\sqrt{Z_{0}}$, illustrating the
relationship between the bond-order variable and the coordination number, see
also Ref. [249]. The bond-order variable $\bar{b}$ typically does not depend
on bond-length, but only on crystal structure. We can therefore read-off the
equilibrium bond length as,

$\frac{a_{0}^{\text{eq}}}{r_{0}}=1-\frac{\ln\bar{b}}{n-m}.$ (102)

Similarly, the energy per bond scales as

$U_{\text{pb}}^{\text{eq}}(\bar{b})=U_{0}\bar{b}^{2-2\mu},$ (103)

explaining the rough structural trends in Fig. 9 for carbon. We can combine
Eqs. (102) and (103) to eliminate $\bar{b}$, yielding

$\frac{a_{0}^{\text{eq}}}{r_{0}}=1-\frac{1}{n}\ln\frac{U_{\text{pb}}^{\text{eq}}}{U_{0}},$
(104)

which is displayed in Fig. 9h.

These analytical equations are useful in parameterization of these potentials.
Albe [297, 298, 300] developed a parameterization procedure, based on data of
ground-state crystalline structures such as the one shown in Fig. 9. The first
step is fitting the parameters $n$ to the universal energy-bond relation given
by Eq. (104). Albe calls the corresponding plot, such as Fig. 9d for copper
and 9h for carbon, a Pauling plot. The remaining (pair) parameters $U_{0}$,
$m$ and $r_{0}$ are then adjusted to the binding energy, bond length and
vibrational frequency of the dimer. In a final step, the parameters that enter
the specific calculation of $\bar{b}$ are fitted to reproduce detailed
structural information in the database.

###  6.2 Defects, melting, and boiling

Reproducing defect energies correctly constitutes an important benchmark for
potentials, as their proper description is required for the atomistic modeling
of thermodynamic and non-elastic mechanical properties. Defects are broadly
characterized into point, line, and planar defects. Reproducing their energies
accurately is more difficult than reproducing an EOS, because atoms near
defects lack the symmetry of ideal crystallographic positions. Parameters
having no effect on the EOS for reasons of symmetry and a small effect on
lattice vibrations can become important for the proper description of defects,
mostly because defects imply a coordination change and might be accompanied by
the generation of radicals or ions in covalently bonded systems.

Vacancies are the most important point defects in simple crystals. Two values
for their energy are associated with it, i.e., the vertical vacancy energy,
$U_{\textrm{vac}}^{(0)}$, for which the lattice in the vicinity of the removed
atom is kept undistorted, and the true defect energy, $U_{\textrm{vac}}$,
which is obtained after lattice relaxation. The $\mathcal{O}(10\%)$ reduction
of the defect energy due to relaxation can be deemed irrelevant when studying
trends, in contrast to the 10 to 105 fold increase in computing and human time
to estimate that correction to within 10%. Different system sizes would have
to be selected, each structure be relaxed to the energy minimum, and results
be extrapolated to the thermodynamic limit, while assessing a vertical defect
energy requires merely one calculation in addition to the ones finding the
crystal’s ground state.

For a short-range potential without explicit angular dependence, the energy
per atom can be estimated with
$U^{\text{eq}}_{\textrm{pa}}=U_{0}Z_{0}^{\mu}/2$ – see also Eq. (99) –, where,
however, $\mu=1$ for pair potentials and $\mu\lesssim 1/2$ for Ducastellium.
Removing an atom from an ideal lattice site and placing it hypothetically into
another ideal lattice site makes its $Z_{0}$ ditched neighbors lose one bond
each, leads to an estimate of,
$U_{\textrm{vac}}=Z_{0}\\{U_{\text{pa}}^{\textrm{eq}}(Z_{0})-U_{\text{pa}}^{\textrm{eq}}(Z_{0}-1)\\}$,
which becomes

$U_{\textrm{vac}}^{(0)}=\frac{U_{0}}{2}Z_{0}\left\\{Z_{0}^{\mu}-(Z_{0}-1)^{\mu}\right\\}\approx\mu
U_{\textrm{pa}}^{\textrm{eq}}$ (105)

where the latter approximation applies if $(1-\mu)/Z_{0}\ll 1$. Thus,
$U_{\textrm{vac}}^{(0)}=U^{\text{eq}}_{\textrm{pa}}$ for a pair potential
($\mu=1$), while for $\mu<1$, our back-of-the envelope calculation provides a
lower, yet good estimate for $U_{\textrm{vac}}^{(0)}$. Taking into account
that the bond lengths of the ditched neighbors is not at the perfect bond
length for their new $Z_{0}-1$ coordination would lead to a relative increase
of $U_{\textrm{vac}}^{(0)}$ of typically less than 5% in a Ducastelle
potential. Fig. 11 confirms the predicted trend that
$\mu^{\prime}=U_{\textrm{vac}}^{(0)}/U_{\textrm{pa}}^{\textrm{eq}}$ is
distinctly less than unity for simple metals, in fact their $\mu^{\prime}$
falls slightly below 0.5, which is the upper bound for $\mu$ approached in the
limit $n\gg m$, mentioned in the text below Eq. (99). Group 14 elements have
substantially larger $\mu^{\prime}$ than metals, however, $\mu^{\prime}$ and
$\mu$ cannot be expected to be similar, whenever bond-angle corrections
matter. Only argon out of the elements considered in Fig. 11 has a vertical
defect energy reasonably consistent with the pair-potential assumption.

Figure 11:  Melting temperature in units of cohesive energy as a function of
$U_{\rm{vac}}^{(0)}/U_{\text{pa}}^{\text{eq}}$ for a variety of elemental
systems condensing in bcc (crosses), fcc (triangles), and diamond-cubic
structure (dc, diamonds). Melting temperatures and cohesive energies were
taken from Kittel’s book [301]. Vertical defect energies are provided for Al
[302], Au [302], Pt [302], Cu [302], Pb [302], Pd [302], Na [303], Ta [302], W
[302], Cr [302]. For C, Si and Ge, $U_{\text{pa}}^{\text{eq}}$ values were
calculated using DFT-LDA. In addition, the vertical defect energy reduction of
Ar was deduced from the ATM parameters from Ref. [149].

One reason why defect energies in metals are important is that they are
believed to correlate with the melting temperature $T_{\textrm{m}}$: Górecki
found the vacancy concentration just below $T_{\textrm{m}}$ in many different
metals to be $n_{\textrm{vac}}^{*}\approx 0.0037$ [304] . Estimating the
relative number of vacancies with $n_{\textrm{vac}}\approx\exp(-\beta
U_{\textrm{vac}})$ would yield
$k_{\textrm{B}}T_{\textrm{m}}/U_{\textrm{vac}}=-1/\ln
n_{\textrm{vac}}^{*}\approx 2/11$, which roughly overestimates the melting
temperatures reported in Fig. 11 by a factor of two. Thus, the simple picture
laid out here is only half the story, or, depending on viewpoint, it is
already half the story. Given the utmost simplicity of the arguments, we would
consider the melting-temperature glass to be half full.

To extend the discussion to the boiling (or sublimation) temperature,
$T_{\textrm{b}}$, the temperature has to be identified at which the Gibb’s
free energy $G$ of the condensed phase equals that of the gas. At low
temperature, $G$ can be crudely estimated with the cohesive energy of the
crystal and at high temperature with $-T_{\textrm{b}}$ times the entropy of an
ideal gas. This gives a quasi-linear relation between
$U_{\textrm{pa}}^{\text{eq}}$ and $T_{\textrm{m}}$ with corrections
logarithmic in $T_{\textrm{m}}$ and $p$. When checking data for metals, we
found them to obey a $k_{\textrm{B}}T_{\textrm{b}}\approx
0.06~{}U_{\textrm{pa}}^{\text{eq}}$, surprisingly well. Unfortunately, we did
not find a name for this relation, although a similar rule certainly exists in
the literature. Even water does not stray too far from this finding with
$k_{\textrm{B}}T_{\textrm{b}}/U_{\textrm{pa}}^{\textrm{eq}}\approx 0.054$,
when expressing its cohesive energy of
$U_{\textrm{pm}}\approx$0.6\text{\,}\mathrm{eV}$$ [305] in ice Ih per molecule
(pm) rather than per atom.

Given the discussions on melting and boiling temperatures, their ratio
$r_{\textrm{mb}}\coloneqq T_{\textrm{m}}/T_{\textrm{b}}$ must be expected to
correlate with the exponent $\mu$. For our representative pair-potential
element argon, melting and boiling temperature differ only by 4% at
atmospheric pressure. In contrast, $r_{\textrm{mb}}(\textrm{alkali})\approx
0.3$, $r_{\textrm{mb}}(\textrm{earth alkali})\approx 0.6$, and most transition
metals being somewhere in between.

Of course, much more detailed analyses are required to estimate reliably
melting temperatures and boiling points from defect and cohesive energies than
the crude arguments presented above. The trouble is, as with potentials, that
improving estimates even marginally requires a substantially increased effort.
Thus, for the sake of brevity, we cannot go into the required level of detail
and instead invite the reader to conduct the fun/instructive exercise of
melting and vaporizing a two-dimensional crystal in a simulation cell having
many times times the crystal’s volume: one time with a simple pair potential
and one time with a generic Ducastelle or related potential. If done
correctly, interesting insights can be gained, for example, on long-range
positional and orientational order in two spatial dimensions, or the
exceedingly small densities required to make the metal gas consist of monomers
rather than of dimers, revealing one more reason why the above presented
reason for the correlation between $T_{\textrm{b}}$ and $U_{\textrm{c}}$
should be bad. For reasons of brevity, we do not only abstain from a more
detailed analysis of the effects of the $U_{\textrm{pa}}^{\text{eq}}\propto
Z^{\mu}$ scaling, but also from reporting trends for surface or grain boundary
energies, although both are central to the properties of materials.

###  6.3 Elastic properties

Deforming a solid from its equilibrium geometry requires the application of
external stress. For a macroscopically affine deformation it can be deduced
from elastic tensors. A brief introduction to them including their relation to
pair potentials is given in Sect. 3.5. While we have already elaborated on the
Cauchy relations, it seems much less appreciated that the pair-potential
approximation also places rather tight bounds on $C_{11}/B$.

The bulk modulus $B$ is easy to evaluate given the equation of state, such as
Eq. (95),

$B=-v_{\text{pa}}\frac{\partial p}{\partial v_{\text{pa}}}.$ (106)

The first constraint to be discussed is the ratio $C_{44}/B$ for fcc crystals
interacting with short-range potentials. Using the results compiled in Eq.
(35), a first estimate of $C_{44}/B=3/4$ is obtained using the Cauchy relation
$C_{44}=C_{12}$ and $B=(C_{11}+2C_{12})/3$ valid for cubic systems interacting
through pair potentials. This is already close for to the value of
$C_{44}/B\approx 0.8$ for fcc Lennard-Jonesium and, by fortuitous error
cancellation, even closer to the experimental value for argon.

There are two leading-order corrections to the nearest-shell approximation,
which both decrease $C_{11}$ relative to $C_{12}$ and thus $C_{44}$ so that
the exact value for $C_{44}/B$ is increased compared to the nearest-shell
approximation. First, the effective spring constant $k_{1}$ associated with
the next-nearest neighbor shell is slightly negative, because next-nearest
neighbors are located well past the inflection point of the LJ potential at
zero external pressure and even further beyond the point at which $k(r)$
becomes negative. Numerically, we obtain $k_{0}\approx 129\,U_{0}/r_{0}^{2}$
and $k_{1}\approx-4.8\,U_{0}/r_{0}^{2}$ when using $a_{0}/r_{0}=0.971$ for fcc
Lennard-Jonesium. This accounts for more than 50% of the difference between
the nearest-shell approximation and the exact value for fcc Lennard-Jonesium.
Second, the dominant effect of more distant neighbors is similar to an
isotropic compressive stress $p_{\textrm{c}}$, pushing the nearest neighbors
more deeply into the repulsion and leading to a reduction of $C_{11}$ by
$p_{\textrm{c}}$ and an increase of $C_{12}$ and thus $C_{44}$ by the same
amount, see also Eq. (39).

The second constraint to be discussed for a short-range potential is the
stability of the bcc phase. In the nearest-neighbor approximation,
$C_{11}=C_{12}$, which means that the tensor of elastic constants is not
positive definite thereby violating the Born mechanical stability criteria
[142]. Including the second shell helps stabilizing bcc, but only if the
interaction is sufficiently long ranged so that $k_{1}$ is positive. For a
Morse potential, we find $C_{12}$ to be less than $C_{11}$, once $n\lesssim 9$
and the binding energy per atom $E_{0}(\textrm{bcc})>E_{0}(\textrm{fcc})$ for
$n\lesssim 6.1$, at which point not even 20% of the interaction energy is
associated with the nearest-neighbor shell. Thus, stabilizing bcc over fcc
using a two-body potential straying not too far away from a physically
meaningful, _atomistic_ potential does not appear to be possible. This can
differ for positively charged colloids repelling each other through Yukawa
potentials [160].

Some solid-state physics textbooks argue that a positive next-nearest neighbor
coupling $k_{1}$ stabilizes the simple cubic structure against shear, polonium
being the single element in the periodic table adopting it as ground state.
However, for $k_{1}$ to be positive at a next-nearest-neighbor distance of
$r_{1}=\sqrt{2}\,a_{0}$, the Morse exponent has to be as small as $n\approx
2.34$, which is much smaller than any reported value for a Morse type
potential that we have come across.

Given the general form of an EAM or TB2M potential of Eq. (57), elastic
constants of simple crystals can be found quite easily with the formalism laid
out in Sects. 3.5,

$\displaystyle v_{\text{pa}}C^{U\eta}_{\alpha\beta\gamma\delta}$
$\displaystyle=$
$\displaystyle\sum_{s}\left\\{2\tilde{U}_{2}^{\prime\prime}(a_{s}^{2})+F^{\prime}(\rho_{0})\tilde{f}^{\prime\prime}(a_{s}^{2})\right\\}a_{s}^{4}\nu_{s}^{\alpha\beta\gamma\delta}$
(107)
$\displaystyle+F^{\prime\prime}(\rho_{0})\sum_{s,t}\tilde{f}^{\prime}(a_{s}^{2})\tilde{f}^{\prime}(a_{t}^{2})a_{s}^{2}a_{t}^{2}\nu_{s}^{\alpha\beta}\nu_{t}^{\gamma\delta},$

where $F^{\prime}(\rho_{0})$ and $F^{\prime\prime}(\rho_{0})$ are the first
and second derivative of $F(\rho)$ evaluated at $\rho_{0}$, respectively,
$\tilde{U}_{2}(S)=U_{2}(r)$, $\tilde{f}(S)=f(r)$, and $S=r^{2}$. In the
nearest-neighbor approximation, this expression simplifies to

$v_{\text{pa}}C^{U\eta}_{\alpha\beta\gamma\delta}=\left\\{2\tilde{U}_{2}^{\prime\prime}(a_{0}^{2})+F^{\prime}(\rho_{0})\tilde{f}^{\prime\prime}(a_{0}^{2})\right\\}a_{0}^{4}\nu_{0}^{\alpha\beta\gamma\delta}+F^{\prime\prime}(\rho_{0})\left(\tilde{f}^{\prime}(a_{0}^{2})a_{0}^{2}\right)^{2}\nu_{0}^{\alpha\beta}\nu_{0}^{\gamma\delta}$
(108)

Similar results for elastic tensor elements have been identified numerous
times [46, 45] including a derivation using the shell-tensor concept [192],
however, mostly for stress-free references. The violation of the Cauchy
relations results exclusively from the last summand on the rhs of Eq. (108),
e.g., for systems with inversion symmetry,

$C^{U\eta}_{12}-C^{U\eta}_{66}\equiv
C^{U\eta}_{1122}-C^{U\eta}_{1212}=\frac{F^{\prime\prime}(\rho_{0})}{v_{\text{pa}}}\left\\{\tilde{f}^{\prime}(a_{0}^{2})a_{0}^{2}\nu_{0}^{11}\right\\}^{2}$
(109)

in the nearest-neighbor approximation. Thus, if the embedding function
$F(\rho)$ has a positive curvature at $\rho_{0}$, as is the case for a
Ducastelle potential with $F(\rho)\propto-\sqrt{\rho}$, then $C_{12}>C_{44}$
follows for cubic systems. The symmetrized shear modulus,
$G\coloneqq(C_{44}+C_{\textrm{s}})/2$, where the so-called tetragonal shear
modulus is defined as $C_{\textrm{s}}\coloneqq(C_{11}-C_{12})/2$, can be shown
to obey $G=\alpha\mu B$ in a nearest-shell Ducastelle potential, where
$\alpha=9/8$ for fcc and $\alpha=1$ for bcc [140]. Thus, the most generic EAM
or TB2M potential predicts $G/B$ to scale linearly with $\mu$, just like the
vacancy energy. In fact, simple metals obey this trend, as demonstrated in
Fig. 12. Nonetheless, this correlation must be taken with a grain of salt,
since bcc is not elastically stable in the nearest-neighbor approximation.

Figure 12:  Vertical defect energy in units of cohesive energy as a function
of the symmetrized, dimensionless shear modulus $G/B$ for a variety of
elemental systems condensing in bcc (crosses), fcc (triangles), and diamond-
cubic structure (dc, diamonds). Elastic constants data for all chosen elements
but Ar were taken from Ref. [68]. Data for $C_{44}/C_{12}$ for Ar from Ref.
[306]. Vertical defect energies are provided for Al [302], Au [302], Ni [302],
Cu [302], Li [303], K [302], V [302], Nb [302]. The vertical defect energy
reduction of Ar was deduced from the ATM parameters from Ref. [149]. For C,
Si, Ge and Sr, vertical defect energies were produced using DFT-LDA, with
$4\times 4\times 4$ super-cell to minimize the interaction between a vacancy
site and its image. The full and dashed gray lines show the prediction of the
Ducastelle potential for fcc and bcc, respectively.

In contrast to the embedding term, angular terms as those used by Keating,
Stillinger-Weber, or Tersoff, counteract shear, so that they lead to a Cauchy
violation having the opposite sign as that in Ducastellium and related
potentials.

###  6.4 Plasticity

Plasticity and ductility hinge on the ability of solids to introduce low-
dimensional structural defects, in particular dislocations, whose presence is
needed to prevent metals from brittle fracture under loading. The more easily
(non-planar) defects are inserted into a crystal, the more ductile it will be.
In a very microscopic picture, metals can be argued to be ductile, because
“delocalized electrons allow metal atoms to slide past one another without
being subjected to strong repulsive forces that would cause other materials to
shatter” [307]. However, going directly from the delocalized nature of the
electrons in metallic bonds is somewhat of a big step and does not provide any
guidelines allowing the degree of ductility to be assessed.

In this section and the literature reviewed therein, vacancy energies in units
of the cohesive energy and shear elastic constants in units of the bulk
modulus are found to be relatively small if interactions are isotropic and
bonds weaken with increasing coordination number. Both ratios turn out to be
of similar order, and, generally less than 1/2 for most metals. Given that the
Poisson’s ratio of an isotropic material satisfies $\nu=(3/2-G/B)/(3+G/B)$, it
can be seen that the low limit of $G/B\to 0$ bring the Poisson’s ratio close
to that of a liquid, i.e., to $\nu=1/2$, while the lower limit for $\nu$ is
just a little above the value valid for pair potentials $\nu=1/4$. Thus, the
smaller $G/B$, the more “liquid-like” the elastic tensor, which, of course,
also implies that the energy to introduce a dislocation, which is linear in
$G$, is relatively small. In fact, the elemental solids with a small
$\mu^{\prime}\mathrel{\mathop{\ordinarycolon}}=G/B$ ratio, but also the ones
with a small $U_{\textrm{vac}}/U_{\textrm{c}}$, are by and large more ductile
than the ones with a larger $\mu^{\prime}$. A nice example is niobium, which
happens to be quite ductile despite condensing in bcc, which, unlike fcc,
lacks close-packed glide planes. Of course, assessing the effect of defects on
the plasticity of a material requires the knowledge of many more properties
than the dimensionless ratios $G/B$, $k_{\textrm{B}}T/U_{\textrm{vac}}$, and
so on [308, 309], however, the trends reflected in our examples, which were
not cherry picked, can certainly be understood from these numbers.

The objective of much atomistic simulation work on metals is to reproduce
thermal and mechanical properties of specific metals and their alloys [35,
309]. The attempt to rationalize trends has become a secondary objective,
which, however, is central to our overview article. Since we are not aware of
a work asking the question how the choice of the potential affects the outcome
of a typical tensile-load experiment given a certain micro structure, we took
the liberty and produced such data. The purpose is to demonstrate that the
trends, which are discussed in regard to ductility, are indeed borne out. To
this end, we melted two-dimensional Ducastellium ($\mu=0.36$, typical for
copper, with $N\approx 33,000$ atoms) and quenched it to small temperature,
whereby crystalline grains of different orientation were produced. The such
produced samples, were then subjected to a tensile loading, where the velocity
of a few hundred atoms in the outermost layer were constrained to assume a
constant value in the direction of loading, and a uniform deformation was
simultaneously applied to the rest of the sample. The exercise was repeated
for Lennard-Jonesium, in which case the configurations were rescaled and
allowed to relax before loading. The orientation of the grains remained
unaltered during the relaxation. Results are shown in Fig. 13.

Figure 13: Comparison of two in-silico tensile tests. Both simulations started
from identical microstructures. The left row assumes Lennard-Jones
interactions, the right - a generic model originally proposed by Ducastelle,
which later formed the basis for embedded-atom and Finnis-Sinclair potentials.
Parameters of Ducastelle potential were taken from Ref. [223]. Simulations
were run using LAMMPS [310, 311]. Colors indicate the orientation of grains,
i.e., the angle $\phi_{i}$ obtained when averaging
$\Phi_{i}\exp(\mathrm{i}6\phi_{i})\mathrel{\mathop{\ordinarycolon}}=\langle\exp(\mathrm{i}6\varphi_{ij})\rangle$
over all atoms within a cutoff, where $\varphi_{ij}$ denotes the angle of a
bond w.r.t. $x$-axis. Atoms, for which, $\Phi_{i}<0.75$, are depicted in
black.

Fig. 13 reveals very clearly that the pair-additive Lennard-Jonesium shows the
characteristics of brittle fracture, while the isotropic, non-pair-additive
Ducastellium produces a force strain relation akin of a ductile metal. In
addition, the pair-additive Lennard-Jonesium breaks along existing grain
boundaries without a change of internal grain boundary structure, while
Ducastellium undergoes massive reconstruction, the precise nature of which
cannot be predicted, as it changes with the random seed for the thermostat. Of
course, these two-dimensional simulations are merely a cartoon for real
solids, among other reasons because stress-induced recrystallization in two
dimensions occurs much more easily than in three, due to a small (hyper-)
surface to volume ratio.

##  7 Charge-transfer potentials

It is sometimes insufficient to treat atoms as if they had a constant
(integer) charge. This calls for methods allowing charge assignments to be
made on the fly. Adjustable coefficients entering the charge-transfer part of
a potential must be gauged on given reference data as any other force-field
parameter. One possibility is to identify them such that inter-molecular
forces or interaction energies from a training set are matched as closely as
possible, which, however, is rarely advisable since errors in certain parts of
the regular part of the potential may be compensated by false charge
assignments, which then jeopardizes the transferability. This is why
reasonable, directly determined reference charges are needed in the training
of charge-transfer potentials.

One reason why knowing partial charges accurately – whatever this means – is
important, is that ionic interactions, so they are present in a system, are on
par with covalent interactions. This is easily seen by comparing their
magnitude. The covalent bond in H2 has a binding energy of 4.74 eV. The
Coulomb interaction in a NaCl molecule with a bond length of $a_{0}=2.36$ Å is
about 6.1 eV, after subtracting sodium’s ionization energy,
$I_{\textrm{Na}}=5.14$ eV, and adding chlorine’s electron affinity,
$A_{\textrm{Cl}}=3.61$ eV, an energy gain of 4.57 eV remains. This needs to be
contrasted to typical LJ interaction parameters of order 0.01 eV. Thus, even
minor changes in partial charges easily alter energies more than relatively
large changes of LJ prefactors.

###  7.1 Determination of reference charges

In the absence of periodic boundary conditions, a charge distribution can be
rigorously characterized in terms of its net charge, dipole moment, quadrupole
tensor, and so on. Nonetheless, there is no unique way to assign _atomic_
charges, _atomic_ dipoles, etc., in a many-body system. Yet, the concept of
formal and partial charge is central in chemical physics and physical
chemistry. As a consequence, many different charge-assignment schemes have
arisen to achieve the ill-defined. In fact, any pragmatist simply _must_
attempt to assign atomic charges allowing interatomic forces to be computed
accurately while providing values that satisfy chemical intuition. For
example, when trying to match the vibrational frequency, the bond length, and
the binding energy of a NaCl molecule with a three-parameter potential, be it
Mie or simple Buckingham with two parameters for repulsion plus one
independent partial charge, the latter turn out in the vicinity of unity so
that chemical intuition is satisfied. Even more, the such obtained parameters
allow quite reasonable predictions to be made, for example, for the lattice
constant, bulk modulus and in fact, all three independent elastic tensor
elements of NaCl in the rocksalt structure, even if truly accurate potentials
require (many-body) corrections. We note in passing that producing numbers
supporting this claim is a great exercise for students and instructors alike.
At the same time, deducing charges from molecular dipoles, which yields
$Q_{\textrm{Na}}=0.75$ e when ignoring the induced dipole on the Cl- ion and
making students deduce that a bond in NaCl to be 25% covalent is a task, which
we deem halfway between useless and insane, even if countless chemistry
lectures and quite a few text books provide discussions along such lines. The
above back-of-the-envelope calculation simply neglects the high polarizability
of anions [90, 93].

Complications in assigning partial charges arise as soon as partial charges
are not close to an integer value, which happens when atoms with
electronegativity differences of, say, $\Delta I\gtrsim 0.5$ eV, are present,
simple alkali halides being the major exception to this rule. Unfortunately,
very few partial charges are assigned as easily as those for sodium or
chlorine atoms in the NaCl molecule. As argued at the beginning of Sect. 7, it
is usually beneficial to separate charge assignment from the construction of
the potential. Existing methods to compute reference charges can be crudely
categorized into those in which atomic charges are deduced from an analysis of
(i) a representation of the wave function (Mulliken [312] and Löwdin [313]),
(ii) the electronic density (Bader [314] and Hischfeld [315]), (iii) the
electrostatic potential (ESP), (iv) displacement-induced changes in
polarization (Born [75]), and frequently forgotten (v) experiment [316]. These
methods will be sketched next. For more detailed overviews, see Refs. [316,
317, 318].

The Mulliken [312] definition of the partial charge starts from the
representation of the total wave function as a linear combination of atomic
orbitals. The charge associated with each orbital is assigned to the atom,
however, results are far from unique, in particular when using large basis
sets. Löwdin [313] removed a certain degree of ambiguity from Mulliken’s
definition that originated from the non-orthogonality of the functions
spanning the basis set. However, uncertainties remain. Bader [314] partitions
the space into regions such that each atom is assigned the charge density in
its vicinity up to the points at which the charge density passes through a
local minimum. However, this procedure ignores the possibility of atomic
charge distributions to interpenetrate. It yields finite charges for the
promolecules, whose charge density is defined to be the superposition of the
atomic charge densities [319]. Hirshfeld [315] originally assumed each atom to
own space with a weight proportional to the electronic density of the neutral
atom, which, however is clearly inappropriate for ions. For example, in an
alkali hydride, H- would be assigned a smaller radius than the metal cation.
Thus, in more advanced methods, the atomic charge densities are replaced with
radially symmetric weighting function that are determined self-consistently
[315, 320]. In such “iterative stockholder approaches”, it remains challenging
to identify good trade-offs between the ESP accuracy and the transferability
of the charges [321]. In conventinal ESP schemes, atomic charges are adjusted
to best match the ESP a save distance away from the nuclei, e.g., outside the
atomic van-der-Waals radii. This becomes problematic for condensed phases, in
which such safe distances are sparse or even non-existent, but also for large
molecules with hidden atoms. Even worse, the ESP appears to be ill-defined up
to a constant in periodic systems so that partial changes can only be gauged
on ESP differences [322]. Born effective charges reflect the change of
polarization occurring in response to atomic displacements, which, for
periodic structures would be analyzed in terms of the Fourier transform of the
dipole moment at small wave numbers. However, since both dipoles and
displacements are vectors, Born charges are not scalars but tensors of rank
two, unless the system is so highly symmetric that all eigenvalues of the
tensor are identical. In addition, Born charges can cast changing atomic
dipoles as effective charge, which can lead to charges exceeding the formal
oxidation number. Central ways to deduce partial charges from experiment [316]
include (i) the matching of electrostatic multipoles, which ultimately is a
limiting case of an ESP method, (ii) Bader-type analysis of the charge density
deduced from x-ray diffraction patterns, but also (iii) a combined analysis of
phonons and dielectric response functions.

Quite a few authors claim to have identified the most meaningful way to assign
partial charges. This can scarcely be true, since partial charges are ill-
defined so that the optimal charge-assignment scheme depends on the relative
weight of target properties. Moreover, the optimum method can also depend on
the symmetry of the problem. Problems arising due to “falsely” assigning an
induced dipole on an atom as a transferred charge vaporize for highly
symmetric structure, in which the lowest-order allowed multipole can be a
hexadecupole. Unfortunately, separating atomic dipoles from charge transfer is
even difficult for diatomic molecules, which could otherwise serve as a well-
defined reference system on which charge-assignment schemes can be gauged.

An optimum choice for a charge to be used in a force-field based simulation
may also depend on whether polarizable dipoles are included or neglected so
that the optimum reference-charge-determination method may depend on the
precise nature of the force field wanting to be parametrized. As a
consequence, we are neither willing nor able to express an opinion if any of
the sketched methods is all in all the optimum choice and believe that it has
to be made case by case.

###  7.2 Regular charge-equilibration approaches

In the ground-state of a full quantum-mechanical system, the functional
derivative of the energy $E$ w.r.t. the electron density $\rho(\mathbf{r})$
disappears, because $\rho(\mathbf{r})$ minimizes $E$. In a coarse-grained
description of the electron density in terms of a set of atomic charges,
$\\{Q\\}$, this minimization principle translates into what Mortier and
varying co-authors [323, 324] called the electronegativity equalization
principle. Assuming that $E$ can be expanded into powers of atomic charges,
the ground state energy can be written as

$\displaystyle U$ $\displaystyle=$ $\displaystyle
U(\mathbf{Q}=0)+\sum_{i}\frac{\partial U}{\partial
Q_{i}}\,Q_{i}+\frac{1}{2}\sum_{i,j}\frac{\partial^{2}U}{\partial Q_{i}\partial
Q_{j}}\,Q_{i}\,Q_{j}+...$ (110a) $\displaystyle=$ $\displaystyle
U(\mathbf{Q}=0)+\sum_{i}\chi_{i}\,Q_{i}+\frac{\kappa_{i}}{2}\,Q_{i}^{2}+\sum_{i,j\neq
i}\frac{J_{ij}}{2}\,Q_{i}\,Q_{j}+...,$ (110b)

while nuclear coordinates remain fixed. Charges adopt themselves such that
they minimize $U$ usually with the constraint of a fixed total charge. In Eq.
(110), the electron density still minimizes the total energy for any
$\\{Q_{i}\\}$, albeit with the constraint imposed through the assigned partial
charges. Thus, the parameters $\kappa_{i}$, $\chi_{i}$, and $J_{ij}$ are
implicitly defined from a mathematical point of view in Eq. (110), even if
their numerical values depend on the method used to relate charge density and
partial charges.

The chemical meaning of the parameters $\kappa_{i}$, $\chi_{i}$, and $J_{ij}$
is best ascertained in limiting cases. Consider atom $i$ to be neutral and
isolated. If one electron is added to it so that $Q_{i}=-e$, the energy
decreases by its electron affinity $A_{i}$, while it increases by its
ionization energy $I_{i}$ when one electron is taken away from it. Thus,

$\displaystyle I_{i}$ $\displaystyle=$
$\displaystyle\frac{\kappa_{i}}{2}\,e^{2}+\chi_{i}\,e$ (111a) $\displaystyle
A_{i}$ $\displaystyle=$ $\displaystyle\frac{\kappa_{i}}{2}\,e^{2}-\chi_{i}\,e$
(111b)

so that the so-called chemical hardness [325] turns out to be
$\kappa_{i}=(I_{i}+A_{i})/e^{2}$ and the electronegativity [326, 327]
$\chi_{i}=(I_{i}-A_{i})/e$. Meanwhile, $J_{ij}$ can be associated with the
Coulomb potential given the charges $Q_{i}$ and $Q_{j}$ of two distant atoms
$i$ and $j$:

$J_{ij}=\frac{1}{4\pi\varepsilon_{0}\,r_{ij}}\textrm{ for }r_{ij}\to\infty.$
(112)

The effect of an external electrostatic potential would have to be added to
$\chi_{i}$. The bold hope now is that partial atomic charges can be predicted
using Eq. (110) and parameterizations, which do not stray too far away from
Eqs. (111) and (112).

Rappé and Goddard [328] pioneered the use of the electronegativity
equalization principle for the on-the-fly determination of atomic charges in
molecular simulation and renamed the approach to charge-equilibration (QEq)
method. To keep charges in (partially) ionic systems from blowing up, they
suggested to shield the Coulomb interactions at short distances through two-
electron Coulomb integrals of the form

$J_{ij}=\frac{1}{4\pi\varepsilon_{0}}\,\int\,\mathrm{d}^{3}r_{i}\,\mathrm{d}^{3}r_{j}\frac{\rho(\mathbf{r}_{i})\,\rho(\mathbf{r}_{j})}{r_{ij}},$
(113)

where the charge densities $\rho(\mathbf{r})$ are those associated with
$s$-type Slater orbitals. This way, the Hessian, i.e., the matrix formed by
$\partial^{2}U/\partial Q_{i}\partial Q_{j}$ remains positive definite even
when atoms approach each other closely. Nonetheless, their claim that QEq
“leads to charges in excellent agreement with experimental dipole moments and
with the atomic charges obtained from the electrostatic potentials of accurate
ab initio calculations” can be seen skeptically: Can Coulomb shielding be
responsible for charges in the alkali metal halides to be so darn close to
unity? And how can QEq reproduce the experimental observation that partial
charges of atoms and ions in the gas phase adopt integer multiples of the
elementary charge? In fact, QEq obviously predicts the ions of a neutral,
dissociated NaCl molecule to each carry a charge of magnitude
$(\chi_{\textrm{Cl}}-\chi_{\textrm{Na}})/(\kappa_{\textrm{Cl}}+\kappa_{\textrm{Na}})\approx
0.4$ elementary charges when truncating the expansion of Eq. (110) after the
quadratic term [329].

QEq has multiple practical shortcomings, which all originate from the ease
with which _partial_ charge can be transferred from one atom to another one
even over large distances. To name a few, it leads to the wrong dissociation
limit of molecules [329], DFT having related issues, QEq makes the
polarizability of polymers grow superlinearly rather than linearly in the
chain length [330, 331], and, similarly, solids adopt the dielectric response
function of a metal [332], i.e., the dielectric constant of a QEq solid is
infinitely large. Finally, QEq makes the dipole of alcohols be linear in the
length of the hydrocarbon chain [333].

One fundamental problem of Eq. (110) is that the curve of lowest average
energy of an _isolated_ system versus the number of electrons turns out to be
a series of straight lines when generalizing the Hohenberg-Kohn theorem to
fractional particle numbers [334]. This property should be reflected in any
coarse-grained approximation to DFT, which is built on the Hohenberg-Kohn
theorem. However, Eq. (110) fails to do so.

Despite the fundamental issues of QEq, metals can be discussed, since partial
charges can be readily transferred over large distances. In a parallel-plate
capacitance geometry, the length over which an external-field induced charge
density decays into the solids scales approximately with $\sqrt{\kappa}$
[332]. In other words, QEq implicitly contains corrections to the behavior of
ideal metals. Thus, mirror charges, or, rather the charge distributions
producing the same electrostatic field outside of the metal as mirror charges
do, are more smeared out in a QEq solid than in an ideal metal. To reflect
this, the term $\kappa_{i}$ to be used in a QEq simulation of metals may best
be parameterized through an appropriate choice of the Thomas-Fermi screening
length, which can benefit, e.g., the description of the interfacial
interactions between a metal and an electrolyte [335]. For nanostructured
electrodes, the finite density of states at the Fermi level leads to an
additional contribution to the capacity of the device [336], which can be
modeled by replacing the quadratic term in Eq. (110) with an appropriate
nonlinear function [337].

###  7.3 Charge-equilibration methods for non-metallic systems

Various ideas have been pursued to suppress the non-local charge transfer in
QEq for dielectric systems. One of the first and most influential approaches,
coined the fluctuating-charge (fluc-Q) model [338], imposes a charge
neutrality constraint on individual molecules. Unfortunately, this bug fix
does not remedy the superlinear scaling of the polarizability with the chain
length of polymers. Moreover, the bonding topology of molecules needs to be
defined making the simulation of bond breaking and formation difficult to
describe.

To suppress long-range charge transfer, Chelli et al. [330] introduced atom-
atom charge transfer (AACT) variables $q_{ij}=-q_{ji}$ so that the charge of
an atom is given by

$Q_{i}=\sum_{j}q_{ji},$ (114)

where $q_{ji}$ is the charge transferred from atom $j$ to atom $i$. In the
AACT method, the bond hardness, $\kappa_{ij}$, replace the atomic hardness
used in regular QEq so that a charge transfer costs a potential energy of
$\kappa_{ij}\,q_{ij}^{2}/2$. The AACT model remedies all shortcomings in QEq
approaches originating from non-local charge transfer, since the AACT method
allows the latter to be (completely) suppressed with (infinitely) large bond
hardnesses. However, for the Hessian of the AACT model to be positive
definite, the bond hardness terms have to be made so large that the dielectric
constants of condensed phases, $\varepsilon_{\textrm{r}}$, turn out to be
barely exceed unity [339]. Moreover, AACT suppresses the observed chain-length
dependent polarizability of short oligomers [331, 332] and, similarly,
produces a zero screening length in dielectrics [332]. Last but not least,
AACT cannot be used to model the polarization of metals that occurs in
response to an electric field.

Reintroducing the atomic hardness term to the AACT model leads to a new model,
which Nistor et al. [339] first meant to call fluc-Q 2. Ultimately, they found
the term split-charge equilibration (SQE) to be more appropriate, since the
model is a split between two models and a “split charge”, $q_{ij}$, formerly
known as AACT variable, is split between two atoms. Thus, their model reads

$U=U(\textrm{QEq})+\sum_{i,j>i}\frac{\kappa_{ij}}{2}\,q_{ij}^{2},$ (115)

where the $Q_{i}$ entering $U(\textrm{QEq})$, i.e., Eq. (110), are defined as
in Eq. (114). Including finite bond hardness suppresses the superlinear
scaling of molecules with linear size [340], as is shown in Fig. 14. The
observable underestimation of the polarizability is supposedly due to the
neglect of atomic polarizability.

Figure 14:  Polarizability of (a) electronegativity equilization model (EEM)
(equivalent to the QEq model) and (b) SQE as a function of the polarizability
deduced from quantum mechanical calculations on a variety of molecules. Green,
red, and blue data points reflect the smallest, the medium, and the largest
eigenvalue of the polarizability tensor for different molecules (including
varying conformations). Reprinted from Verstraelen, T., Van Speybroeck, V.,
and Waroquier, M. J. Chem. Phys. 131, 044127, (2009), with the permission of
AIP Publishing.

Other charge-transfer models exist ensuring neutral separation limits. In one
line of approach, chemical potential differences between distant atoms are
screened with functions motivated from overlap integrals between orbitals on
different atoms [329]. However, this procedure implicitly claims that the
altitude difference of two objects becomes less when moving them apart
laterally, which is not meaningful. In contrast, a very appealing approach is
a model in which electrons are treated as shells, which are not tight to one
individual atom [156], whereby, in principle, history dependence can be
mimicked along the lines of the NaCl dissociation described in the
introduction. There are also approaches using non-local softness matrices, the
arguably most advanced one being the atom-condensed Kohn-Sham approximated to
second order (ACKS2) [341]. Its purpose is twofold: First, defining the
expressions for the softness matrix, which is the inverse of the hardness
matrix, from expressions accessible from DFT calculations, and second to
employ the result in simulations as a mean to deduce partial charges.
Ultimately, ACKS2 contains SQE model as a limiting case. The downside of ACKS2
is that a $2N\times 2N$ matrix must be defined, where $N$ is the number of
charges, while the number of split charges is only of order $ZN$. A frequently
raised argument against SQE is that the large bond hardnesses occurring at
large separation cause numerical difficulties. However, they can be addressed
using appropriate pre-conditioners in square-gradient minimization or large
split-charge inertia in extended Lagrangian approaches.

One of the main limitations of the original SQE model is that it cannot
produce isolated ions as the $\kappa_{ij}$ diverges at large $r_{ij}$. This
limitation can be overcome by introducing oxidation numbers $n_{i}$ so that an
atomic charge is given by

$Q_{i}=n_{i}\,e+\sum_{j}q_{ji}$ (116)

with $n_{i}\in\mathbb{Z}$ in a so-called redox-SQE formalism [79, 159]. The
remaining part of the SQE formalism is kept, except that force-field
parameters, including chemical hardness and electronegativity, must now be
assigned individually for each oxidation number. The redox-SQE approach allows
a system to effectively move on different potential energy surfaces (if you’re
a chemist) or different Landau-Zener levels (if you’re a physicist). Jumping
between them requires the oxidation number of one atom to be increased by an
integer, typically by one, while that of another, nearby atom must be
decreased by the same number. The corresponding discrete dynamics can be
implemented in practice, for example, in terms of Monte Carlo dynamics.

###  7.4 Properties resulting from charge-transfer potentials

A central aspect of polarizable and charge-transfer potentials is to reproduce
the dielectric response of (condensed) media correctly, the most important
quantity from a continuum perspective being the dielectric constant
$\varepsilon_{\textrm{r}}$. The relation between point-dipole polarizability
and $\varepsilon_{\textrm{r}}$ is text-book material and reflected in the
Clausius-Mossotti (CM) relation, Eq. (42), while that of a bond-polarizable
model without point dipole polarizability was worked out by Nistor and Müser
[332]. Combining the treatment of point and charge-transfer dipoles does not
appear to have been presented in the literature hitherto. However, it can be
easily achieved when following Hannay’s derivation of the CM relation [342].
It starts from the insight that the full electrostatic field of a dipole
contains a contribution
$\mathbf{E}_{\textrm{int}}=-\mathbf{p}\,\delta(\mathbf{r})/(3\varepsilon_{0})$
acting inside of it in addition to its usually stated far field.

Since each charge leads to a mean electrostatic field of zero, the mean
electrostatic field of a charge distribution is also zero, however, outside
the dipoles, it is reduced by $\mathbf{E}_{\textrm{int}}$ averaged over the
volume that each point dipole occupies. Thus, the average field that a test
charge sees outside of dipoles is the external electric field
$\mathbf{E}_{\textrm{ext}}$ (which is assumed to be constant) minus
$\langle\mathbf{E}_{\textrm{int}}\rangle=-\rho\,\mathbf{p}/(3\varepsilon_{0})$,
where $\rho$ is the density of point dipoles. For a rocksalt lattice, or, when
discretizing a homoegeneous material into cubes of linear size $a_{0}$ and
arranging the local charge transfer variables, or, split charges in a vector
$\mathbf{q}$, the dipoles obey associated with the charge-transfer,
$a_{0}\mathbf{q}$ and the point dipoles obey

$\displaystyle\frac{a_{0}\mathbf{q}}{\kappa_{\textrm{s}}}$ $\displaystyle=$
$\displaystyle\mathbf{E}_{\textrm{ext}}+\frac{\rho}{3\varepsilon_{0}}\,\mathbf{p}$
(117a) $\displaystyle\alpha\,\mathbf{p}$ $\displaystyle=$
$\displaystyle\mathbf{E}_{\textrm{ext}}+\frac{\rho}{3\varepsilon_{0}}\,\mathbf{p},$
(117b)

respectively. The $x$-component of $\mathbf{q}$ is the split charge donated
from an atom or grid point to its right neighbor, and so on, while
$\kappa_{\textrm{s}}$ is the corresponding bond hardness. The total
polarizability is nothing but the total dipole density
$\mathbf{P}=\rho(\mathbf{p}+a_{0}\,\mathbf{q})$, while
$\varepsilon_{\textrm{r}}$ is defined through
$\mathbf{P}=(\varepsilon_{\textrm{r}}-1)\varepsilon_{0}\mathbf{E}_{\textrm{loc}}$
within linear response. After some minor algebra

$\varepsilon_{\textrm{r}}-1=\frac{\rho/(a_{0}\,\kappa_{\textrm{s}})+\rho\alpha}{\varepsilon_{0}-\rho\alpha/3}$
(118)

can be obtained, which contains the CM relation in the limit of
$\kappa_{\textrm{s}}\to\infty$ as well as
$\varepsilon_{\textrm{r}}-1=\rho/(a\kappa_{\textrm{s}})$ valid for bond-
polarizable models ignoring point dipoles [332]. Once bond polarizability is
ignored in a charge-transfer model, $\varepsilon_{\textrm{r}}$ diverges, which
means that the system effectively turns metallic.

Continuum corrections to the above treatment would require the wave-number
dependence of the charge-dipole and dipole-dipole coupling but also the
chemical hardness coupling to be taken into account. The latter would augment
the bond hardness in the presented treatment with $\kappa(qa_{0})^{2}$. A
wave-number dependent dielectric constant and thus a dielectric smearing or
decay length $\delta$ ensues. For a pure bond-polarizable model and the
geometry considered above,
$\delta=\sqrt{\varepsilon_{0}a_{0}\kappa/(1+\varepsilon_{0}a_{0}\kappa_{\textrm{s}})}$
[339]. One of the consequences of this result is that ideal mirror charges do
not accurately represent the dielectric response of real metals [335].

The potential importance of modeling dielectric response correctly is also
demonstrated in Fig. 15. It shows two pairs of two dislike solids, or,
clusters, which can exchange charge, specifically, oxidation numbers, when
they are sufficiently close to each other. In reality, this distance would be
a weak function of the approach and retraction velocity. In one pair, both
materials are metals, in the other, they are insulators. A lot of “electrons”
transfer in the metallic couple and swapped charge delocalizes over the entire
samples, mostly on their surfaces. In the insulating couple, only one charge
exchanges, which remains localized near the former contact points after
retraction.

Figure 15: Contact electrification with history-dependent potentials between
two metals (left) and two dielectrics (right). Reprinted by permission from
Springer: Eur. Phys. J. 86, 337, Towards time-dependent, non-equilibrium
charge-transfer force fields, Dapp, W. B. and Müser M. H., [COPYRIGHT] (2013)

Note that Fig. 15 is a proof-of-principle demonstrator for the redox-SQE
method [79]. It does not allow one to conclude that contact electrification
between dielectrics primarily occurs through electron transfer. There are many
good reasons to believe that ion exchange matters most [343].

##  8 Machine-learned potentials

In the past two decades, machine learning has conquered the world of
interatomic potentials, a development that is ongoing and accelerating. The
main idea is to equip a machine learning (ML) technique with a large database
of structures for which energies, forces, stresses, etc. are computed using
well-converged first-principle methods333See
https://libatoms.github.io/GAP/data.html for links to high-quality _ab-initio_
databases for C, Si, and other materials.. Energies (or sometimes forces) are
then extrapolated (or “predicted”) based on (local) structural similarity.
Structural similarity is either encoded by directly comparing to
configurations stored in this database, or it can be implicit in the weights
of a neural network or the points of a sparse Gaussian process. Recent methods
supplement such computational approaches with experimental data [344].
Additionally, predictions are not limited to energies and forces. Examples of
other properties obtained from ML regression are core energies and X-ray
photoelectron spectra [345, 346, 347, 348].

Machine-learned potentials (MLPs) are sometimes characterized as being
parameter-free, however, one may also claim that the number of parameters is
gigantic. Occam’s razor suggests that models should have few parameters. Yet,
overfitting of high-dimensional models can be avoided by using
parameterization procedures that are based on sound statistical principles,
such as Bayesian inference. Additionally, all physically-justified potential
forms (with few parameters) are based on a chain of approximations, e.g.,
going from the Born-Oppenheimer approximation, to DFT, neglecting three-body
integrals in the tight-binding approximation, neglecting nonorthogonality and
so on in the bottom-up construction of bond-order potentials for
semiconductors. It is often unclear to what extend the final form is still a
good approximation for interatomic interactions. To say it with the words of
Bazant, Kaxiras and Justo [239]: “ _In the case of Si, the abundance of
potentials in the literature illustrates the difficulty of the problem and
lack of specific theoretical guidance. In spite of the wide range of
functional forms and fitting strategies, all proposed models possess
comparable and insufficient overall accuracy. It has proved almost impossible
to attribute the successes or failures of a potential to specific features of
its functional form._ ”

MLPs target reproducing the high-fidelity calculations used to construct the
database. In that sense, they are higher accuracy (with respect to the
reference DFT database) than more empirical or semiempirical models, but this
may mean a lack of transferability to structures not explicitly considered in
the database. This pitfall is of course not limited to MLPs: As an example,
the first incarnation of the Tersoff potential for silicon [243] turned out to
have bcc and not diamond cubic as the ground state [349]. Yet, potentials
constructed on carefully curated databases give excellent transferability even
to obscure crystal phases, see e.g. Ref. [350] for a current example, again
for elemental silicon. This highlights the specific importance of database
quality for MLPs.

Most MLPs consist of three basic components: First, a database of structures,
already outlined above. Second, a “descriptor” that encodes the chemical
environment of an atom into a unique fingerprint. In the absence of external
fields, this fingerprint must be objective, i.e., it must be invariant with
respect to translation and rotation. It must also be invariant with respect to
the exchange of atoms of the same species. This allows the data set to be much
reduced compared to the situation, where the ML technique has to learn
objectivity by itself. An additional complication is that the database itself
may not be perfectly objective, due to numerical errors in the underlying
electronic structure calculations. Since we want forces, fingerprints also
must be continuously differentiable at least once, and ideally, they also have
continuous second derivatives. The third component is a regression framework,
the “learning” aspect of the potential. Two broad directions of development
have occurred here: Using neural networks, which can be regarded as a model
having many parameters, or, Gaussian processes, which can be regarded as a
scheme directly extrapolating between points in the database. The regression
framework takes the fingerprints as input and predicts per-atom energies as
output.

Before continuing with more details, we would like to point out that many
reviews by researchers more qualified than us to discuss MLPs have appeared
recently, see Refs. [351, 29, 31, 32, 352]. This review focuses on simple
functional forms that allow analytical calculations to be made. Yet, not
discussing at least the basic ideas behind MLPs could be a severe omission in
present times. More importantly, some of the discussion highlighting
similarities between classical and MLPs might help the development of either
one.

###  8.1 Chemical fingerprints

Chemical fingerprints come in many forms, and we here discuss just a select
few. The interested reader is referred to other literature, e.g. Refs. [353,
352], for more comprehensive descriptions of the wide variety of fingerprints.

The arguably first and one of the simplest example of a chemical fingerprint
consists of the coordination number and the bond-conjugation variable employed
by Brenner in his potential for hydrocarbons [10, 251, 11]. Given a bond
$i$-$j$, Brenner used the coordination numbers $N_{i}$, $N_{j}$ and a bond
conjugation parameter $N_{ij}^{\text{conj}}$ as a fingerprint for the state of
the bond (see Sect. 5.7 for more details). These variables are made continuous
by computing the coordination numbers with a smooth cutoff,
$N_{i}=\sum_{j}f_{\text{c}}(r_{ij})$ where $f_{\text{c}}(r)$ is one of the
cutoff functions discussed in Sect. 9.1. This fingerprint is then used to
correct the bond-energy of a Tersoff potential through a lookup table. The
lookup procedure continuously interpolates between integer values of the
triplet $(N_{i},N_{j},N_{ij}^{\text{conj}})$ using splines. The table entries
were fitted to a small database of experimental properties of hydrocarbon
molecules.

Brenner’s tables are used to correct a baseline potential, a procedure that is
sometimes called a $\Delta$-potential, e.g. in Ref. [32]. Other descriptors
are needed for full-fledged machine-learned energies. As a simple instructive
example, we consider two and three-body terms in the sense of Eq. (2). An
objective descriptor for two-body interaction is simply the bond distance
$r_{ij}$. For a three-body term, we can use the triplet of distances
$(r_{ij},r_{ik},r_{jk})$ as a descriptor. While this descriptor satisfies
objectivity, it still has permutational symmetries with respect to the
exchange of the positions of like atoms $k$ and $j$. Bartók and Csányi [351]
suggested the use of the descriptor

$\mathbf{q}_{ijk}=(r_{ik}+r_{ij},(r_{ik}-r_{ij})^{2},r_{jk})$ (119)

which is equivalent to $(r_{ij},r_{ik},r_{jk})$ but symmetric with respect to
the permutation $j\leftrightarrow k$.

We note that common potentials, like Axilrod-Teller-Muto (Sect. 4.1) and
Stillinger-Weber (Sect. 5.4), which in their original formulation depend on
the $(r_{ij},r_{ik},r_{jk})$ triplet, can also be written in terms of the
descriptor $\mathbf{q}_{ijk}$. These analytical potentials obey the above
permutational symmetries by construction. However, in a machine-learning
context, where an energy is essentially obtained by extrapolating from a given
set of database entries, there is no guarantee that such permutational
symmetry will be fulfilled by the extrapolated quantities. This requires
reformulation in a form such as Eq. (119).

The specific descriptor of Eq. (119) was used in the Gaussian approximation
potential (GAP) [354] for amorphous carbon [355], later extended for improved
accuracy of crystalline phases [356]. This potential followed a specific
three-step construction: In the first step, a (nonparametric, i.e. machine
learned) pair-potential was fitted to a database using the scalar $r_{ij}$ as
a descriptor. The second step involved correcting this potential using a
three-body term employing the descriptor given by Eq. (119). Finally, this
potential was again corrected using the many-body SOAP descriptor described
further below. Incorporating simple descriptors (including a pair descriptor)
into the construction of the potential gave improved transferability. This
basic potential without the SOAP descriptor can be regarded as a machine-
learned variant of the cluster expansion, see Eq. (2), in an embedding medium.

A successful strategy for obtaining more generic descriptors is the use of
atom-centered symmetry functions introduced by Behler and Parrinello [357,
358]. The atom-centered fingerprints consist of radial symmetry functions

$q_{i}^{(1,\alpha)}=\sum_{j\not=i}e^{-\eta^{(\alpha)}(r_{ij}-R_{s}^{(\alpha)})^{2}}f_{\text{c}}^{(\alpha)}(r_{ij}).$
(120)

These functions are Gaussians centered at a pair distance $r_{s}^{(\alpha)}$,
which are forced smoothly to zero at a cutoff distance via
$f^{(\alpha)}_{\text{c}}(r)$. They are similar to coordination numbers, but
put specific weight on atoms at the distance $R_{s}^{(\alpha)}$. In addition
to these radial functions, there are angular functions

$\displaystyle q_{i}^{(2,\beta)}$ $\displaystyle=$ $\displaystyle
2^{1-\zeta^{(\beta)}}\sum_{j,k\not=i}(1+\lambda^{(\beta)}\cos\vartheta_{ijk})^{\zeta^{(\beta)}}$
(121) $\displaystyle
e^{-\eta^{(\beta)}(r_{ij}^{2}+r_{ik}^{2}+r_{jk}^{2})}f_{\text{c}}^{(\beta)}(r_{ij})f_{\text{c}}^{(\beta)}(r_{ik})f_{\text{c}}^{(\beta)}(r_{jk}).$

We have added the superscripts $(\alpha)$ and $(\beta)$ to emphasize, that
those functions are not evaluated for a single parameter vector
$(\eta,R_{s},\lambda,\zeta)$ and cutoff function but for a set of these
parameters. The fingerprint vector $\mathbf{q}_{i}$ is then constructed by
combinations of $q_{i}^{(1,\alpha)}$ and $q_{i}^{(2,\beta)}$ over these
parameters sets. This allows one to construct fingerprint vectors of arbitrary
dimension but those fingerprints may contain redundant information. In the
following, we will refer to the components of the fingerprint vector as
$q_{i}^{(\mu)}$ and drop the double notation $(1,\alpha)$, $(2,\beta)$,
implying that the index $\mu$ also refers to the type of symmetry function and
not just its parameters.

To emphasize connection to classical potentials, the $q_{i}^{(1)}$ look like
the onsite $\xi_{ii}$ terms occurring in our general functional form, Eq.
(54), while the $q_{i}^{(2)}$ look like offsite $\xi_{ij}$ terms. Indeed, if
we replace the Gaussian in Eq. (120) with an exponential, the expression
becomes identical to the local density of Ducastellium, see Eq. (59b).
Similarly, $q_{i}^{(2)}$ are constructed like the many-body contributions in
the SW and MEAM potentials.

Another important strategy for obtaining chemical fingerprint is the smooth
overlap of atomic positions (SOAP) descriptor, introduced by Bartók, Kondor
and Csányi [359]. It can be regarded as a generalization of the bond-
orientational order parameter of Steinhardt, Nelson and Ronchetti [360]. The
idea is to map the atomic configurations onto an atomic density field
$\rho(\mathbf{r})$, which is given by atom-centered functions,

$\rho(\mathbf{r})=\sum_{i}\delta(\mathbf{r}-\mathbf{r}_{i}),$ (122)

where $\delta$ is broadened to a Gaussian with a characteristic length
$\ell_{\text{at}}$. Note that $\rho(\mathbf{r})$ is invariant to permutations
of atoms by construction and resembles the density constructed for EAM
potentials. To obtain angular information, this density field is then expanded
into a set of spherical harmonics $Y_{lm}(\theta,\phi)$ for each atom $i$,

$\rho_{i}(\mathbf{r})=\sum_{nlm}C_{i,nlm}g_{n}(r)Y_{lm}(\theta_{i},\phi_{i})$
(123)

with local angles $\theta_{i}$ and $\phi_{i}$ of the vector
$\mathbf{r}-\mathbf{r}_{i}$. The functions $g_{n}(r)$ are orthogonal radial
basis functions. The power-spectrum of the angular coefficients $C_{i,nlm}$,

$p_{i,nn^{\prime}l}\propto\sum_{m}C_{i,nlm}^{*}C_{i,n^{\prime}lm}$ (124)

is objective and hence fulfills symmetry considerations. Selected elements of
this power spectrum, typically chosen up to an upper value of $n$, form the
SOAP descriptor.

SOAP taken to higher body order – and with $\ell_{\text{at}}\to 0$ as in Eq.
(122) – forms the basis of the atomic cluster expansion (ACE) [361, 362]. We
recommend study of Ref. [361] also because it outlines the relationship
between current machine-learning approaches and more traditional EAM and bond-
order potentials. First parameterizations of (performant) ACE potentials are
starting to appear in the literature [363].

###  8.2 Regression

Once we have a fingerprint vector, we need to interpolate (or since we are in
a high-dimensional space, rather extrapolate) atomic energies. This is the
task of the regression framework. Regression frameworks take as input a
fingerprint vector and predict _per-atom_ energies $U_{i}$, such that the
total energy is simply $U=\sum_{i}U_{i}$. This decomposition corresponds to
the one used in all other potentials discussed in this review.

Regression can be as simple as a linear (ridge) regression. This is for
example employed in spectral neighbor analysis potential [364], which uses the
fingerprints introduced in the original GAP [354]. A more flexible regression
strategy is the use of Gaussian process regression [253], which lends its name
to the GAP potential [354]. An advantage of Gaussian processes is that there
is a clear statistical interpretation of Gaussian process regression, which
allows one to predict confidence intervals on energies in addition to
expectation value or maximum probability estimates. In an active learning
framework, this can help to identify regions in a database for which
additional first-principles calculations (or experiments) are necessary [365].
The final expression for the per atom energies takes the simple form,

$U_{i}=\sum_{t}\alpha_{t}K(\mathbf{q}_{i},\mathbf{q}_{t}),$ (125)

where the sum over $t$ runs over all elements of the training database. The
fingerprints $\mathbf{q}_{i}$ of the current environment of atom $i$ are
compared to all fingerprints $\mathbf{q}_{t}$ in the training database. The
function $K(\mathbf{q}_{i},\mathbf{q}_{t})$ is the kernel that encodes a
measure of distance in fingerprint space (and implicitly includes a mapping
into a high-dimensional representation, typically called “feature space” in
the ML literature [253]). From Eq. (125) it becomes immediately clear that the
numerical effort of a GAP is linear in the number of points in the reference
database. It is therefore crucial to _sparsify_ the database, i.e., to remove
entries that are close together in fingerprint space and do not improve the
predictive power of the potential. Sparsification is highly nontrivial and we
refer to the original literature on details [32].

A key feature of Gaussian process regression is that a notation of smoothness
of the potential energy landscape is build-in. A common choice for the kernel
function $K$ is a squared exponential (Gaussian) [355, 253], which containsa
set of distance scales $\ell_{\alpha}$ as a hyper parameter ,

$K(\mathbf{q}_{i},\mathbf{q}_{t})\propto\exp\left\\{\frac{(q_{i,\alpha}-q_{t,\alpha})^{2}}{2\ell_{\alpha}^{2}}\right\\}.$
(126)

The final potential energy is smoothed over $\ell_{\alpha}$, which avoids
parasitic local energy minima or artifacted phonon spectra, which are
sometimes found when regressing with neural networks [366], or even when
simply interpolating with splines [222]. GAPs additionally make use of a dot-
product kernel for the many-body (SOAP) contribution to the energy [355]. The
dot-product kernel has no intrinsic length-scale, but the SOAP descriptors
contain intrinsic smoothing through the width $\ell_{\text{at}}$ of the
underlying atom-centered Gaussians. The expansion coefficients $\alpha_{t}$
are obtained using a Bayesian regression scheme that is described in detail in
the pertinent literature [253].

An interesting use of Gaussian process regression is the on-the-fly
parameterization [367] of a force field based on electronic structure
calculations carried out on snapshots throughout a molecular dynamics
trajectory. Li, Kermode and de Vita [368] used a small database of reference
configurations, generated from past simulation steps of a molecular dynamics
trajectory, to extrapolate forces to the future. Their scheme differs
qualitatively from the ones discussed above by predicting forces rather than
energies. The reason for this is that unlike the per-atom energy, a per-atom
force is an observable that can be computed straightforwardly within first-
principles techniques for a subset of atoms in a system. This leads to a
nonconservative force field, yet use of this method to accelerate a quantum-
mechanical region in concurrent multi-scale modeling [367], e.g., of fracture
[369, 370], introduces no additional complication as the coupling scheme
itself is typically non-conservative.

We note that GAPs are linear models (in feature space), as the final energies
are linear combinations of Gaussians. ACE potentials [361, 371] introduce an
additional nonlinear mapping. Using a square-root (as in EAM or TB2M
potentials) reduced the number of basis functions required for high-accuracy
regression for copper [371]. Fully nonlinear regression can be achieved with
neural networks. Behler and Parrinello [357] employed the specific form

$U_{i}=f_{\text{a}}^{(2)}\left\\{w_{01}^{(2)}+\sum_{\nu}w_{\nu
1}^{(2)}f_{\text{a}}^{(1)}\left(w_{0\nu}^{(1)}+\sum_{\mu}w_{\mu\nu}^{(1)}q_{i}^{(\mu)}\right)\right\\},$
(127)

for regression in their neural network potential. Here,
$f_{\text{a}}^{(1)}(x)$ and $f_{\text{a}}^{(2)}(x)$ are sigmoidal “activation
functions” and $w_{\nu 1}^{(1)}$ and $w_{\nu 1}^{(2)}$ are weights. The
weights are numbers that are adjusted in the process of fitting the neural
network. The sum over $\mu$ runs over all entries of the feature vector
$\mathbf{q}_{i}$ of atom $i$. This network has a single “hidden layer” with
activation function $f_{\text{a}}^{(1)}(x)$ and the index $\nu$ runs over the
outputs of this layer. The original fit for silicon by Behler and Parrinello
employed three nodes in the hidden layer and $48$ entries in the fingerprint
vector, amounting to a total of $196$ weights. Equation (127) resembles the
many-body contribution $U_{\text{m}}$ of classical potentials for open-shell
systems, see Eq. (53).

##  9 Limiting the interaction range

Chemistry is local for the most part [372]. As a consequence, forces between
distant atoms are minor, unless they are charged or carry dipoles. By
neglecting interactions with distant atoms, much computing time can be saved
at the expense of minor systematic errors, which can often be made marginal
with mean-field or other coarse-graining corrections [373, 374, 375].
Unfortunately, the locality of chemistry cannot necessarily be specified in
terms of cut-off radii, because to what extent two atoms “see” each other
depends on various factors. Screening and cutting are somewhat related but the
bonding topology and density must be known to express screening through
cutting.

A metal atom will scarcely be exposed to the charge density of another atom if
a third atom is sandwiched in between them. Baskes incorporated such screening
conceptually into MEAM potentials [226, 376] motivated by the observation that
“ _the second-nearest neighbors do not play a part even in the bcc structure_
,” despite the ratio of nearest to next nearest-neighbor shell being as small
as $1.15$. He continued that “ _this surprising result gives some support to
the procedure we have used here in ignoring all but first-neighbor
interactions_ ” [226] (in the description of crystals). Defects and disorder,
both central to materials properties, place large demands on the assessment of
bonding topology and make it necessary to define gray shades between bonded
and non-bonded. In realistic descriptions, those cannot be judged using only
the distance between two atoms.

Dispersive interactions between noble-gas atoms are not screened so that
limiting their interaction range (without accounting for long-range
corrections) is merely an exercise in finding a good compromise between
accuracy and computational efficiency. Therefore, the task of cutting
potentials might appear technical or even trivial at first sight. Doubling the
cutoff radius $r_{\textrm{c}}$ beyond which interactions are ignored requires
roughly an eightfold computational effort in three spatial dimensions to
construct neighbor lists and to evaluate pair potentials, which can be
unjustifiable if the total dispersive interaction is already accurate.
However, undesired qualitative artifacts can arise when cutting improperly.
Brittle materials can turn ductile [377, 378] and thereby make it on the cover
of contemporary text books [379], or, nano-tube-based water pumps work
endlessly without external driving due to improper cut-offs [380], the report
of which [381] having made it through the exquisite reviewing process of
Nature journals. Nonetheless, cutting potentials is a prerequisite for the
linear scaling of numerical effort per time step with the number of atoms $N$
and for the efficient parallelization on parallel, high-performance computing
systems using domain-decomposition [382, 383, 384, 385, 386, 310, 311]. Even
long-range Coulomb potentials can be computed with a numerical effort scaling
linearly or quasi-linearly with $N$ [387, 388, 389], which requires atom-atom
distance evaluations to be cut off at a finite distance.

Unfortunately, choosing optimum cut-off radii and procedures depends not only
on the potential but also on the property of interest. Someone studying phase
equilibria between different thermodynamic phases wants the potential depths
to be reproduced as accurately as possible, which can be achieved with a
brutal truncation of the potential at a given $r_{\textrm{c}}$ without
shifting. Implementing such a potential is readily done for Monte Carlo,
however, for a molecular dynamics simulation, one would have to account for
the $\delta$-function singularity in the interaction force at
$r_{\textrm{c}}$, which no sympletic integrator in the world would be pleased
about. On the other hand, someone interested in mechanical properties can care
less about the depth of energy minima but should instead be careful about
potential curvature and the smoothness of the force at the cutoff. A bond
whose tensile force disappears only linearly at $r_{\textrm{c}}$ cannot be
broken adiabatically, which risks to exaggerate loss moduli or viscosities: a
spring attached to a bond would snatch into or out of bonding whenever the
bond length $r$ passes through $r_{\textrm{c}}$, whereby energy is dissipated
[390].

Three more general points are worth discussing before cutting to the chase:
First, when two phases like fcc and hcp have a similar thermodynamic
potential, the choice of cutting can affect what phase is preferred [391].
However, at that point it might always be wise to consider other effects like
many-body interactions and thermal or quantum fluctuations of the nuclei.
Second, changing the cut-off procedure or cut-off radius may require a new
parameterization of the potential. Third, and perhaps most importantly, not
only potentials can be cut or smoothed but all operations discussed here below
can also be applied to forces. However, cutting forces generally leads to non-
conservative forces unless they originate from central, two-body potentials.

###  9.1 Cutting

####  9.1.1 Pair functions

Not only pair potentials must be cut off but also expressions that many-body
potentials depend on, such as the charge density in EAM potentials, or the
interaction between point multipoles. Cutting procedures for pair functions
can be broadly categorized into cut-and-shift and smoothing.

Cut-and-shift procedures [392] redefine a pair potential $U_{2}(r)$ such that
all its derivatives up to $l$’th order approach zero continuously at
$r_{\textrm{c}}$ and remain zero for $r>r_{\textrm{c}}$:

$U_{2,l}^{\textrm{shift}}(r)=\left\\{U_{2}(r)-\mathcal{T}[U_{2}(r),r_{\textrm{c}},l]\right\\}\Theta\left(r_{\textrm{c}}-r\right).$
(128)

Here, $\mathcal{T}[U_{2}(r),r_{\textrm{c}},l]$ is the $l$’th order Taylor
expansion of $U_{2}(r)$ about $r_{\textrm{c}}$.

The effect of cut-and-shift procedures on $U_{2}(r)$ and on the equation of
state is shown exemplarily for a Morse potential ($n=6$,
$r_{\textrm{c}}/r_{0}=2$, $l=0,...,2$) in Fig. 16a. It becomes evident that
“large” $l$ lead to a poor representation of the binding energy unless
$r_{\textrm{c}}$ is large, which, however costs computing time. However,
serious artifacts can occur when potentials are truncated but also when they
are merely shifted with $l=0$. The EOS of an athermal crystal is discontinuous
even for $l=0$ as revealed in Fig. 16b. Thermal fluctuations smear out and
potentially eliminate this discontinuity. Nonetheless, caution remains
advisable. In principle, a $l$’th-order shifted potential can have a
discontinuity in $\partial^{l}p/\partial V^{l}$ so that a force-shifted
potential risks to cause a discontinuity in the compressibility or any elastic
tensor element including its imaginary, i.e., viscous part, as just discussed.

Figure 16: Left panel: Truncated (black), cut-and-shift (blue), shifted-force
(red), and shifted-curvature (green) Morse potential with exponent $n=8$ and a
cut-off distance of $r_{\textrm{c}}=2.0\,r_{0}$. The cut-and-shift procedures
include $l=0$ (blue), $l=1$ (red), and $l=2$ (green). Right panel: energy per
atom of an ideal fcc crystal as a function of the nearest-neighbor distance
$a_{0}/r_{0}$ using the shifted potential (blue), shifted force (red), and
shifted curvature (green) for the Morse potential with $n=8$ and a cutoff
radius $r_{\textrm{c}}=2.0\,r_{0}$. Dots and squares denote reduced volumes at
which a neighbor-shell radius is equal to $r_{\textrm{c}}$.

An alternative to cut-and-shift potentials is the use of cut-off or smoothing
functions $f_{\text{c}}(r)$, which are multiplied with the true pair potential
to yield

$U_{2}^{\textrm{cut}}(r)=U_{2}(r)f_{\text{c}}(r).$ (129)

The general property of smoothing functions is that they are unity below an
inner radius, $r_{\textrm{i}}$, which may be zero or even negative. They
decrease past $r_{\textrm{i}}$, approach zero at $r_{\textrm{c}}$, and remain
zero ever after. Formally they can be expressed as

$f_{\textrm{c}}(r)=\Theta(r_{\textrm{i}}-r)+\varphi_{\textrm{c}}(x)\Theta(r-r_{\textrm{i}})\Theta(r_{\textrm{c}}-r)$
(130)

with $x=(r-r_{\textrm{i}})/(r_{\textrm{c}}-r_{\textrm{i}})$. Many of the
lessons learned for cut-and-shift applies to smoothing, in particular w.r.t.
what artifacts are produced depending on the order $l$ with which derivatives
of the potential disappears at $r_{\textrm{c}}$. However, an additional risk
of smoothing functions is that using “large” $l$ does not guarantee artifacts
to be suppressed in particular when $r_{\textrm{i}}$ is close to
$r_{\textrm{c}}$. Prominent smoothing functions $\varphi_{\textrm{c}}(x)$ are:

$\displaystyle\left\\{1+\cos(\pi x)\right\\}/2$ $\displaystyle
l_{\textrm{i}}=1;\;l_{\textrm{c}}=1$ Tersoff [243] (131a)
$\displaystyle(1-x)^{3}(1+3x+6x^{2})$ $\displaystyle
l_{\textrm{i}}=2;\;l_{\textrm{c}}=2$ MEAM [225] (131b)
$\displaystyle\left\\{8+9\cos(\pi x)-\cos(3\pi x)\right\\}/16$ $\displaystyle
l_{\textrm{i}}=2;\;l_{\textrm{c}}=2$ Murty & Atwater [393] (131c)
$\displaystyle\exp\left\\{-(2x)^{m}\right\\}$ $\displaystyle
l_{\textrm{i}}=m-1$ (unknown) (131d)
$\displaystyle\exp(\lambda+\lambda/(x-1))$ $\displaystyle
l_{\textrm{i}}=0;\;l_{\textrm{c}}=\infty$ Stillinger & Weber [228] (131e)
$\displaystyle 2(1-x)^{m}/(1+(1-x)^{n})$ $\displaystyle
l_{\textrm{i}}=0;\;l_{\textrm{c}}=m-1$ Mishin [217, 218], (131f)

where integer Mishin exponents are commonly selected as $m=n=4$. In Eq. (131),
$l_{\textrm{i}}$ and $l_{\textrm{c}}$ specify which highest $l$’th-order
derivative of the smoothed function approaches zero at $r_{\textrm{i}}$ and
$r_{\textrm{c}}$, respectively. Moreover, $l_{\textrm{i}}=0$ means that
$\varphi_{\textrm{c}}(r_{\textrm{i}})=1$ while
$\varphi^{\prime}_{\textrm{c}}(0)\neq 0$, while omitting $l_{\textrm{c}}$
means that $\varphi_{\textrm{c}}(x=1^{-})\neq 0$. Recently, it was found that
cutting with polynomials constructed such that
$l_{\textrm{i}}=l_{\textrm{c}}+1$ can simultaneously improve binding energies
while reducing cut-off artifacts [394]. This is because cut-off induced
discontinuities of (the derivatives of) short-range potentials are distinctly
larger at $r_{\textrm{i}}$ than at $r_{\textrm{c}}$.

####  9.1.2 Cutting Coulomb potentials

Wolf et al. [395] suggested that the Coulomb interaction in an overall neutral
system can be cut off without producing systematic errors. The idea is to
represent the Coulomb interaction through the Ewald summation [396], keep the
$q=0$ contribution (which could be seen as a kind of mean-field correction),
and to ignore any contribution at non-zero wave vector $\mathbf{q}$. Some
authors find support for the correctness of Wolf’s optimistic assessment [286,
397], while others don’t [398, 399]. In fact, small cutoff radii can be used
for homogeneous systems when using appropriate cut-off functions, e.g.,
clearly less than $10$ Å in silica melts [394]. However, problems arise when
structural heterogeneities extend on length scales $\lambda>r_{\textrm{c}}$
[394]. They induce undulations in chemical potentials with wave vectors of
order $q=2\pi/\lambda$, which favor a charge separation. Coulomb interactions
suppress that separation with terms scaling as $1/q^{2}$. The pertinent
restoring forces are neglected for the most part in the Wolf summation.

Another argument for why the Wolf summation cannot be cutoff in heterogeneous
systems can be made when considering a gold nugget separated by a large
distance from a chunk of lithium. If electrons can equilibrate, gold will
acquire a negative charge due to its larger electronegativity. While the
charge density on each lump will decrease with their size, the absolute force
attracting them will increase. This mechanism also matters for discharge
simulations of batteries or related devices [159], despite the presence of a
Helmholtz double layer near the electrodes. If a Maxwellian demon were to keep
all electrolyte atoms in place and closed an electrical switch allowing some
charge to flow between anode and cathode, a current would flow, which would
extend much beyond an elementary charge. Thus, local charge neutrality cannot
be obeyed during the entire time of the discharge process. Once the demon
releases the electrolyte particles, discharge can continue after the
electrolyte rejuvenated the double layers. Modeling the described process with
a Wolf summation properly would require one to use a cut-off exceeding the
anode-cathode separation. Thus, the Wolf summations, or better generalizations
thereof avoiding a discontinuous force at $r_{\textrm{c}}$ [286, 394], can be
made for homogeneous systems, but steer clear of it otherwise, in particular
in any biophysical context [287].

####  9.1.3 Many-body potentials

Smoothing functions are not only applied to pair potentials but also in many-
body potentials, sometimes with the goal to imitate screening. In a
crystalline structure, interactions with atoms in the nearest-neighbor shell
should be unmodified but those with the next-nearest neighbors should be
screened – often those are approximated to be fully screened. Figure 17
illustrates the effect of cutting the interaction in a Ducastelle potential at
various distances. When cutting with a distance-based criterion, energies
(Fig. 17a) are forced to drop to zero somewhere between first and second
neighbor shell of the equilibrium crystal (since we do not want interaction
with second neighbors at equilibrium - they are screened). This leads to a
force or pressure overshoot during homogeneous volumetric dissociation of the
crystal, illustrated in Fig. 17b. Attempting to break such a crystal with an
external stress overestimates the critical stress for fracture by a factor of
four. For this reason, some works have adopted cutting forces rather than
energies for modeling fracture [400], despite the fact that this leads to a
nonconservative force-field. The proper solution to this issue is the use of
screening functions, discussed in the next section. They yield the true
nearest-neighbor potential (shown by the black line in Fig. 17) and hence
smooth dissociation of the crystal.

Figure 17: (a) Energy per atom as a function of distance and (b) pressure as a
function of volume for a Cu fcc crystal described by a Ducastelle potential
[223]. The MEAM cutoff function given by Eq. (131b) was applied to the
repulsive potential and the density. $U_{\rm{c}}$ and $a_{0}$ are the cohesive
energy and the nearest-neighbor distance for an fcc crystal described by a
Ducastelle potential, respectively.

###  9.2 Screening

In order to reflect screening through screening rather than through cutting,
Baskes et al. [226, 376] proposed to replace cutoff functions with
environment-dependent screening functions, $S_{ij}$, which themselves are
products of three-body functions $0<S_{ijk}<1$,

$S_{ij}=\prod_{k\not=i,j}S_{ijk}.$ (132)

The $S_{ijk}$ are constructed such that $S_{ijk}=0$ when atom $k$ screens the
interaction between $i$ and $j$ completely, while $S_{ijk}=1$ when atom $k$
does not screen the $ij$ interaction at all. Thus, the screening function acts
like a topological cutoff function, $0\leq S_{ij}\leq 1$. Baskes proposed
different functional forms for $S_{ijk}$. The refined expression is [376]

$S_{ijk}=\begin{cases}1&\text{if}\quad C_{ijk}\leq C_{\text{min}}\\\
\exp\left\\{-\left[(C_{\text{max}}-C_{ijk})/(C_{ijk}-C_{\text{min}})\right]^{2}\right\\}&\text{if}\quad
C_{\text{min}}<C_{ijk}<C_{\text{max}}\\\ 0&\text{if}\quad C_{ijk}\geq
C_{\text{max}}\end{cases},$ (133)

which drops from unity to zero between $C_{\text{max}}$ and $C_{\text{min}}$.
The quantity $C_{ijk}$ characterizes the geometry of the ellipse that passes
through the three atoms $i$-$j$-$k$. It is given by

$C_{ijk}=\frac{2\left(X_{ik}+X_{jk}\right)-\left(X_{ik}-X_{jk}\right)^{2}-1}{1-\left(X_{ik}-X_{jk}\right)^{2}},$
(134)

the square of the ratio of the lengths of the two half-axes of that ellipse.
Here $X_{ik}=r_{ik}/r_{ij}$ is the distance between atoms $i$-$k$ normalized
by the length of the central $i$-$j$ bond.

Figure 18 illustrates this concept. The bond is unscreened ($S=1$) if there is
no third atom in its vicinity (Fig. 18a). As a third atom moves into the
shaded ellipsis (Fig. 18b – the geometry of the ellipsis is described by the
values of $C_{\text{min}}$ and $C_{\text{max}}$) the bond becomes screened
($S=0$, Fig. 18c). It becomes clear, that if we simply rescale all coordinates
(as in the homogeneous, volumetric dissociation of a crystal), $S$ remains
constant for each bond and an energy would follow the true nearest neighbor
curve, see the black line in Fig. 17a. The screening function $S$ is hence a
scale-invariant formulation of a “cutoff” procedure.

Figure 18: Illustration of the screening concept. The black atom moves along
the hypothetical trajectory indicated by the shaded atoms and the red arrow.