Interatomic potentials: achievements and challenges

ABSTRACT

Interatomic potentials approximate the potential energy of atoms as a function of their coordinates. Their main application is the effective simulation of many-atom systems. Here, we review empirical interatomic potentials designed to reproduce elastic properties, defect energies, bond breaking, bond formation, and even redox reactions. We discuss popular two-body potentials, embedded-atom models for metals, bond-order potentials for covalently bonded systems, polarizable potentials including charge-transfer approaches for ionic systems and quantum-Drude oscillator models mimicking higher-order and many-body dispersion. Particular emphasis is laid on the question what constraints ensue from the functional form of a potential, e.g., in what way Cauchy relations for elastic tensor elements can be violated and what this entails for the ratio of defect and cohesive energies, or why the ratio of boiling to melting temperature tends to be large for potentials describing metals but small for short-ranged pair potentials. The review is meant to be pedagogical rather than encyclopedic. This is why we highlight potentials with functional forms sufficiently simple to remain amenable to analytical treatments. Our main objective is to provide a stimulus for how existing approaches can be advanced or meaningfully combined to extent the scope of simulations based on empirical potentials.

Graphical abstract

KEYWORDS:

• Many-body potentials
• embedded-atom method
• tight-binding
• bond-order potentials
• charge transfer
• dispersive interactions
• machine learning

1. Introduction

Interatomic potentials are functions of nuclear coordinates approximating the electronic ground state energy, or for metals, the electronic free energy of a system. Forces on atoms, as needed in molecular-dynamics simulations, can be obtained by calculating the gradient of these functions with respect to the nuclear coordinates. Frequently, the terms interatomic potential and force field are used synonymously. There is, however, a subtle and sometimes important difference between the two. Many-body force fields are not necessarily gradients of a scalar function, in contrast to “real” interatomic forces when electrons are given enough time to equilibrate. In this review, we focus on force fields that can be represented as gradients of scalar functions approximating the exact interatomic potentials. Before getting started on details, let us take a step back.

In his Lectures on Physics, Richard Feynman wondered what statement would contain the most information in the fewest words, if all of scientific knowledge were to be destroyed in some cataclysm, and only one sentence could be passed to the next generations of scientists [1]: I believe it is the atomic hypothesis that all things are made of atoms – little particles that move around in perpetual motion, attracting each other when they are a little distance apart, but repelling upon being squeezed into one another. Two-body potentials reflect this generic behavior, except those describing two ions carrying a charge of identical sign, since they repel each other even at large separation.

In fact, quite a bit can be learned from studying two-body potentials. For example, the crystalline structure of many elemental or binary crystals can be rationalized, as it is done in any better text book on solid state physics. Moreover, the functional form of constitutive laws is often independent of the precise nature of the potentials, such as Hooke’s law, which states restoring forces of solids in equilibrium to be linear in small deformations. Even many non-linear constitutive equations do not depend on the details of the potential. The exponents with which compressibility or specific heat diverge as the temperature or pressure approaches the fluid-vapour critical point are identical for all substances irrespective of their specific interactions [2–4]. The way how the viscosity of polymers increases with molecular weight can also be described with the help of two-body potentials, as long as they prevent two polymers from crossing each other [5]. These are but a few examples for the success of two-body potentials. The only thing that does depend on chemical detail, it almost seems, are boring prefactors. Unfortunately, or, depending on your viewpoint, fortunately, this is not quite right.

Realistic parameterizations of two-body potentials, including Morse [6], Buckingham [7], or Lennard-Jones (LJ) [8], favor close-packed assemblies of atoms, e.g. face centered cubic or hexagonal closed packed lattices. Thus, neither molecular crystals as those formed by oxygen or nitrogen at small temperature, nor layered crystals like graphite at ambient conditions could be thermodynamically stable. Even those elements that do like to close pack may not be describable by two-body potentials. For example, the proper description of the elastic properties of metals and their ductility hinges on many-body terms [9]. The importance of directed interactions and thus of many-body terms is particularly apparent for carbon, where depending on the hybridization of individual carbon atoms, different interatomic forces ensue. In the worst case, history dependence of interatomic, or rather, interionic forces can arise. To illustrate this point, assume a NaCl molecule dissociates into ions in water and the water is evaporated later so that two isolated ions emerge. If, however, the NaCl molecule had been separated slowly in an inert gas environment before the nuclear degrees of freedom had been brought to their final destinations, two neutral atoms would have formed. This is because the ionization energy of sodium exceeds the electron affinity of chlorine. Real powerful potentials should mimic such a process correctly and account for charge transfer in the appropriate way, though the notion that forces arise as an unambigous function of nuclear coordinates would have to be abondened. From this discussion, we can see that no chemistry and none of its implication (batteries, laptops, soccer, or life in a more general sense) could arise if atomic systems could be accurately described with two-body potentials. In brief, without many-body interactions, the world in general, and science in particular, would be much less exciting than it actually is.

The root of many-body potentials is that two atoms change their interaction when additional atoms are present, because their electronic structure changes. For instance, two hydrogen atoms prefer to form a strong covalent bond between them rather than to remain lonely. However, as soon as an oxygen enters the scene, the hydrogens stop being homo and happily form a heteronuclear water molecule with the oxygen.

Many-body potentials were much advanced over the last few decades. Realistic, large-scale simulations of several thousand atoms can nowadays be conducted even on commodity computers, at least for selected compounds such as simple metal alloys or hydrocarbons. Force-field based simulations of increasingly complex systems and processes become possible, including those during which atoms rehybridize in the course of a chemical reaction [10–12]. However, the need for further improvement remains, in particular for systems in which two or even more different bonding types (ionic, metallic, covalent, dispersive) combine. This text is meant to provide a rudimentary understanding for why specific functional forms for potentials were chosen and what a given (class of) potential may or may not achieve. With this, we hope to stimulate insight on how to combine or to generalize different potential classes so that systems held together by different bonding types can be better simulated in the future. Since an incredible number of papers on potentials has been published, we are certain to have missed important contributions, even if we did our best to identify the original key literature.

Focusing on concepts leaves us little room to provide readers with the best possible parameterization for a given substance and application. Performance evaluations of potentials and their parameterization containing information like the specific parameterization by author A using potential B reproduces properties C and D of the material E but does a poor job on properties F and G are probably better looked up in one of several excellent and important databases [13–18]. In addition, we wish to refer to many excellent reviews [19–33] or even text books [34,35] on interatomic potentials emphasizing specific parameterization of potentials more than this review. Here, we discuss materials- or element-specific numbers only when we explore the transferability of a potential, which is its ability to accurately predict properties, structures, and/or stoichiometries, which it was not adjusted for. In addition, one system is highlighted for each class of bonding type. To this end, we chose argon as the representative of systems bonding through van-der-Waals, because its dispersive interactions are half way between those of water and CH₂ or CH₃ units in hydrocarbon chains. Carbon, copper, and rocksalt (NaCl) complement the list to represent covalently bonded materials, metals, and ionic systems.

2. Classification, construction, and consequences of (many-body) potentials

2.1. Brief classification of interaction potentials

The way in which interactions are incorporated into potentials is so diverse that they can be classified according to several main criteria: the chemical nature of the bond, whether a bonding topology is prescribed or allowed to change, two-body versus many-body interactions, and the degree of empiricism or level of theory used to construct the functional form of the potential and to gauge adjustable parameters.

Atoms with open valence shells form covalent or metallic bonds, while atoms with closed valence shell interact through van der Waals forces. The latter include repulsion at short separation and dispersive forces, which result from the mutual induction of dipoles and other multipoles owing to quantum-mechanical fluctuations of the valence shell. Each bonding type has its own characteristics how interatomic forces change with interatomic separation and as a function of their environment, which motivates the distinction between open-shell and closed-shell potentials. Van der Waals or “non-bonded” interactions have energies on the order of thermal energy at room temperature and are thus weaker than covalent and metallic bonds. In addition, there can be Coulombic interactions due to partial or close-to-integer charges on atoms or ions. Interactions between atoms can also be a combination of all of the above, the most prominent example involving the hydrogen bond.

Potentials with a fixed bonding topology range from simple bead-spring models [5] to highly sophisticated, chemically realistic valence force fields containing explicit bond-angle or torsional terms [36–43]. Two atoms interact differently with each other, depending on whether they are considered bonded or non-bonded to each other, even if all involved atoms carry the same chemical symbol, or, if they correspond to the same coarse-grained entities, such as a so-called united atom reflecting, for example, a CH₂ or CH₃ group [39]. The design of valence force fields is in a rather mature state [36–43], which is why we will discuss them only peripherally. This does not prevent us from highly recommending their use for any situation, in which bond breaking and bond formation can be ruled out, simply because they are computationally lean while being sufficiently accurate for many purposes. Nonetheless, weaknesses remain, in particular the frequently poor treatment of higher-order and many-body dispersion as well as of charge transfer, which, however, are all best explained assuming free atoms or ions rather than bonded entities as reference. In a biophysical context, a proper treatment of dispersive interactions can be particularly important when molecules transfer between aqueous and hydrophobic environments, since the dispersive interactions with carbon-hydrogen chains are stronger than with water molecules [44]. Adsorbing higher-order or many-body effects by tweaking pertinent Lennard-Jones interaction coefficients is then even more hazardous than for homogeneous phases. Thus, readers being interested in improving bonded potentials for biomolecular simulations have no excuse to stop reading here.

The distinction between two-body and many-body potentials should be self-explanatory. The latter can be further subdivided into explicit or implicit. An explicit many-body potential does not require on-the-fly adjustments of parameters like atomic charges or dipoles. In contrast, implicit many-body potentials necessitate the determination of such quantities. This is commonly achieved by minimizing expressions for the potential energy with respect to degrees of freedom representing a coarse-grained version of the full electronic response. Also explicit many-body potentials can be categorized into further subgroups. They can be cast in terms of analytical functions, or numerical tables, or they can be extrapolated from a large set of high-dimensional data as is done in machine-learned potentials [29,31–33].

Potential classes are sometimes also named differently depending on the quantum-mechanical framework from which they are motivated. For example, embedded-atom potentials [9,45] are best motivated from density-functional theory, while second-moment tight-binding potentials [46,47] arise, as their name says, quite naturally from the tight-binding approximation. Potentials lacking a direct theoretical justification are also called empirical. While we could elaborate much more on how to further classify potentials, we feel compelled to cut to the chase.

2.2. Constructing (many-body) potentials

Formally, it is possible to expand the potential energy of a (classical) many-atom system into two-body, three-body, and higher-order contributions, where represents the positions of atoms, i.e.

(1)

Here, the index in the ‘s on the r.h.s. of the equation gives the order of the many-body expression. For non-elemental systems, the also depend explicitly on the atomic indices, or more precisely on the nature of the atoms interacting with each other. The underlying assumption of the expansion is that electrons are in a well-defined state, e.g. in their ground state or in thermal equilibrium.

As long as external fields are present, two-body and higher-order interactions generally do not obey Galilean or rotational invariance, which is illustrated in Figure 1. It depicts two atoms with full valence shells, which are polarized by the electrostatic field of an external dipole. Depending on the relative position of the external dipole to the two atoms, the atoms will experience an electrostatic repulsion, as in panel (a), or attraction, as in panel (b), in addition to their previous interaction.

Figure 1. Effect of an external dipole, marked by a positive and negative sign, on secondary, induced dipoles indicated by arrows. Depending on the location and orientation of the central dipole, be it permanent or dynamic, the induced dipoles repel (a) or attract (b) each other.

In the absence of external fields, the two-body potential can be reduced to depend only on the interatomic distances , in which case Eq. (1)(1) reduces to

(2)

Expansions of the form given by Eq. (2)(2) are sometimes called cluster potentials [19,35]. Unfortunately, even this Galilean invariant expansion Eq. (1)(1) is of limited (direct) use for systems of practical interest, because its convergence can be extremely slow. This is best seen when considering an ion in close proximity to a metallic surface. The ion induces a charge distribution in the metal that depends on the shape of the metal surface turning a formal expansion into explicit many-body terms useless. Likewise, comparing local structural motives to an existing database, as is done with machine-learned potentials, does not appear to be promising in this example, although they frequently allow to effectively truncate an expansion similar to that outlined in Eq. (2(2) ), however, with atoms not being in vacuum but surrounded by other atoms.

Equation (2)(2) converges relatively quickly only for closed-shell atoms, meaning atoms whose valence shell is filled, like noble gas atoms and singly charged alkali cations or halogen anions [48]. Even then convergence might not be truly impressive. Two-body interactions commonly used to simulate condensed phases of noble gas atoms are effective pair potentials rather than true pair potentials accurately reproducing the second virial coefficient [49]. The difference between effective and true potentials compensates to some extent for omitted many-body effects. Thus, pair and three-body potentials trained exclusively in the bulk risk to be inaccurate when applied to surfaces and gases. Despite its frequent ineffectiveness, Eq. (2)(2) allows the difference between two opposite philosophies to the construction of potential function to be discussed: the bottom-up and the top-down design, which would lead to the construction of non-empirical and empirical potentials, respectively.

In the bottom-up approach, pen-on-paper quantum chemistry, density-functional theory (DFT), or any related electronic structure calculations provides analytical results and/or small-scale data, such as forces on individual atoms as a function of the atomic configurations or energies of diatomic molecules as a function of bond length. Such information can motivate the functional form of a two-body potential or be used to construct soulless look-up tables. From a philosophical point of view, it could be argued that inputting experimental information on dimers, trimers, etc., falls into the bottom-up approach, as the electronic structure problem is solved by nature rather than by computers. Of course, the downside of letting nature do that job is that it is quite difficult to get enough experimental data on small molecules and clusters that would allow the higher-order terms in the series of Eq. (2)(2) to be determined to any significant accuracy. Only pair potentials can probably be deduced to high precision from measurements of the vibrational and rotational spectra without having to make serious restrictions on their functional form.

In the top-down design, no atomic-scale information is provided, but instead collective, typically macroscopic properties, such as elastic properties, equations of state (EOS), surface energies, or the temperature dependence of the specific heat including heats of melting and evaporation. In the sense of an inverse problem, or reverse coarse-graining, potential energy functions are constructed such that the available information is reproduced. Historically, the first attempts to design potentials followed this approach. The arguably most systematic top-down design makes use of the virial expansion [50,51], in which the pressure of a many-particle system in thermal equilibrium and in the absence of external fields is expanded into a power series of the number density ,

(3)

where is the thermal energy, while and denote the second and third virial coefficient, respectively. The so-called cluster expansion [52,53] allows the various virial coefficients to be related to two-body, three-body, and higher-order interactions. For example, knowledge of the second virial coefficient

(4)

allows the pair potential to be reconstructed by an inverse procedure, though a numerically stable inversion works best if is represented as a function with few adjustable parameters. The third virial coefficient, , depends on two- and three-body interactions, and so on. Thus, in principle, two substances may have similar second but different third virial coefficients so that knowledge of allows deviations from pair additivity to be quantified and adjustable parameters of a three-body potential to be gauged.

Unfortunately, the virial expansion does not converge for every state point, due to the existence of thermodynamic discontinuities a.k.a. phase transformations [54]. This is one reason why the expansion of Eq. (2)(2) cannot be fully parametrized from a cluster expansion after all. Yet, early calculations performed by van der Waals in the spirit of a cluster expansion, revealed that the leading-order corrections to the ideal-gas EOS requires the two-particle interaction in simple gases to obey [55]

(5)

at large distances. Today, is called the dipole-dipole dispersion coefficient.

About 50 years after van der Waals, London [56] rationalized why atoms and molecules with closed valence shells attract each other through pairwise additive interactions by coupling the quantum-mechanical ground state fluctuations of their dipoles in the far-field approximation. Extending London’s treatment to higher-order electrostatic multipoles and beyond second-order perturbation theory leads to refinements of the theory, which are outlined in Sect. 4. Despite the existence of corrections, the bottom-up approach reveals quite clearly that the exponent in Eq. (5)(5) is essentially exact, at least as long as electrostatic interactions can be treated as instantaneous (or the speed of light as infinitely large) [57,58]. Tweaking the exponent to better match reference data would quickly result in overfitting of the potential: a better match of the reference data would deteriorate the description of interatomic forces between two isolated particles at large distances. Thus, bottom-up and the top-down design of interatomic potentials are complementary to each other but should converge to similar results assuming sufficient and sufficiently accurate input.

2.3. Consequences of many-body interactions

A common way to assess the relevance of many-body interactions is to determine what percentage of the cohesive energy is related to the exact or to an effective two-body interaction. Such an analysis may be misleading, because defect energies or elastic properties may not be described well even if the total energy of a crystal or the radial distribution function of a liquid reproduces ab-initio or experimental data. We can always adjust a pair potential that reproduces the cohesive energy of a specific crystalline phase (discussed extensively below) or that fits a specific liquid pair distribution function (see e.g. Refs. [59,60]). These properties often benefit from the annihilation of or from the insensitivity to many-body contributions, but the resulting potentials are then not transferable between crystal structures, temperatures, or other states of the system. Even three-body correlation functions may be poor properties on which to gauge potentials. Liquid copper and liquid argon just above their crystallization temperature both assemble in structures similar to that of random-sphere packing [61], despite their potentials being utterly different.

Induced dipoles or charge transfer decrease the total energy (and hence make the “bonds” stronger), or they would not occur. In contrast, most other many-body interactions weaken bonds, i.e. the energy per bond in metals decreases with increasing coordination number , typically with for metals and even more quickly for covalently bonded systems (see Sect. 5 or Refs. [46,62,63]). This scaling has important consequences not only on what crystal structure is assumed but also on the relative increase of boiling relative to the melting temperature, which is roughly 4% for argon but 100% for copper.

Another frequent consequence of many-body terms is that bond lengths become shorter with decreasing coordination or increasing bond order. One of the best known example may be the bond length between two carbon atoms: Å (, ethane, C₂H₆), 1.34 Å (, ethylene, C₂H₄) 1.20 Å (, acetylene, C₂H₂), where states how many hydrogen atoms each carbon is bonded to. Similarly, the spacing between layers near a free metal surface contract [64], while common short-range pair potentials would predict them to expand.

Historically, the non-additivity of atomic potentials lead to the dismissal of the atomic hypothesis among many prominent scientists in the 19ʹth century, see Ref. [65] for a well written and very enlightening account of the debate. The assumption was that atoms, so they exist, should only interact through central potentials and be pairwise additive. Cauchy and Poisson [65] and also potentially St. Vernan [66] had shown that this reduced the number of independent elastic tensor elements, e.g. from 21 to 15 for a triclinic crystal and from three to two for a cubic crystal. Specifically, they found, see Sect. 3.5 for a derivation, that any permutation of the four indices of an elastic tensor element should leave it unchanged so that relations, now known as Cauchy relations, like, for example, , or, (in cubic systems ) in Voigt or Nye notation, would hold. It then came as a blow to the atomic hypothesis when Voigt noticed that none of the crystals he had investigated satisfied the Cauchy relation within experimental errors [67]. Figure 2 epitomizes results on the violation of the pair-potential assumption. The latter is also frequently stated in terms of the Cauchy pressure , which vanishes for athermal, inversion-symmetric, classical crystals interacting with pair potentials.

Figure 2. Ratio of elastic tensor elements and as a function of in units of the bulk modulus for a variety of cubic systems. ( in cubic systems.) The blue line holds for pair potentials, while the red line reflects elastically isotropic materials. The thin gray line assumes the generic model proposed by Ducastelle, introduced in Sect. 5.2, which has one dimensionless parameter affecting the and ratios in the nearest-neighbor approximation. Centrosymmetric lattices (closed circles) include face-centered cubic (fcc), body-centered cubic (bcc), rocksalt (rs), and caesium chloride (cc), while diamond cubic (dc) and zinc blende (zb) lack inversion symmetry (open diamonds). Experimental and simulated data are included for selected crystals (Cu [68], Au [68], K [68], V [68], C [68], Ge [68], NaCl [68], LiD [69], LiF [68], MgS(rs) [70], MgS(zb) [70], ZnTe [71], CsCl [72], SiO₂ [73]) as well as calculated elastic constants for CuZr [74].

Rather than abandoning the assumption of the pairwise additivity of potentials, many scientists rejected the atomic hypothesis alltogether. However, Rutherford’s scattering experiments in the early 20ʹth century removed all legitimate doubt about it, even if quite a few muddleheads keep taking issue with it up to this day in happy concert with deniers of human-made climate change. Finally, Born [75] reconciled the properties of the elastic tensor of real materials with the atomic hypothesis by showing that one (of several) reasons for the breakdown of Cauchy relations is the existence of many-body interactions.

It should certainly not be concluded that solids obeying the Cauchy relations automatically satisfy the pair-potential assumption, since there can be fortuitous or symmetry-induced cancellation of many-body effects. For example, the dipole polarizability of anions in the rocksalt structure cannot reveal itself in the elastic tensor for symmetry reasons.

3. Two-body potentials

In this section, we introduce the most generic two-body potentials for pairs of atoms forming covalent and metallic bonds in the condensed phase as well as potentials for ionic and van der Waals interactions. We discuss these two classes separately, because atoms with open electron shells bond primarily through the sharing of electrons while closed-shell atoms interact predominantly through van der Waals or ionic forces. Both classes reflect Feynman’s mantra that [atoms attract] each other when they are a little distance apart, but repel upon being squeezed into one another [1]. Interatomic potentials with well-motivated functional forms can reproduce the equation of state, reasonably well if their small number of adjustable parameters are gauged on a few accurate reference values. However, parameters are only transferable to other structures or properties if the pair-potential assumption is justified.

An element specific analysis of the potentials will be made for carbon (C), copper (Cu), sodium chloride (NaCl), and argon (Ar), which are representative of covalent, metallic, ionic, and van-der-Waals bonding, respectively. To set the stage for further discussion, their pair-potential energies are shown in Figure 3 as a function of the molecular bond length . The reference data shown therein is not necessarily the most accurate on the market, but their overall errors should be minor compared to those of lean pair potentials, i.e. the data should be sufficiently accurate to test if a given functional form of the two-body potential is justified.

Figure 3. (a) Two-body potentials of different diatomic molecules, , normalized to the binding energy as a function of the bond length in units of the equilibrium bond length . Symbols reflect high-accuracy reference data, while lines show fits of pair potentials introduced in Sect. 3 to the energy minimum. Reference data include the Aziz potential for Ar₂ [76], and quantum chemical calculations for C₂ [77], Cu₂ [215], and NaCl [339]. Numbers used for the fits are listed in Table 1. (b) Two-body potentials with a curvature in the minimum corresponding to that of the Lennard-Jones potential.

Table 1. Comparison of pair-potential parameters obtained after fitting adjustable coefficients to accurate reference data on diatomic molecules (mol) and on crystals, where dc and rs are the diamond cubic and rock salt structure, respectively. Numbers indexed with a star are functions of the adjustable parameters. Used potentials are, Born-Mayer/Buckingham for argon and sodium chloride, and the Morse potential for carbon and copper. Plots of the respective potentials are shown in Figure 3. The first seven neighbor shells are included into the fit.

		(eV)	(Å)		(Å)
Ar	(mol)	0.01234	3.757	14.6	0.26*
Ar	(fcc)	0.01155	3.811	14.1	0.27*
C	(mol)	6.318	1.250	5.82	0.21*
C	(dc)	2.868	1.804	5.34	0.34*
Cu	(mol)	1.823	2.275	6.80	0.33*
Cu	(fcc)	0.3964	2.717	7.81	0.35*
NaCl	(mol)	5.152*	2.442	7.91	0.31*
NaCl	(rs)	7.199*	2.309	7.47	0.31*

3.1. Morse potential

One of the early, great successes of quantum mechanics was Heitler and London’s quantitative description of the chemical bond in a hydrogen molecule in 1927 [78]. Their theory is text-book material and will only be sketched here. Heitler and London linearly combined the atomic wave functions of the two hydrogen atoms forming a H₂ molecule, antisymmetrized the spin-wave function to reflect the Fermi principle, and evaluated the binding energy using first-order perturbation theory, thereby providing a lower bound for the molecular binding energy . Sugiura [79] succeeded in identifying the correct solution for the exchange-energy integrals. Errors of order 40% for and bond length remain, see Figure 4. Treating in the electronic ground-state wave functions as a variational parameter, errors can be reduced by a factor close to two.

Figure 4. Binding energy assuming different quantum chemical approaches to the binding in H₂. Full lines are fits to that data near the minimum assuming the functional form proposed by Morse. Dashed lines represent cubic fits to the minimum. The inset shows the same data but normalized to the respective equilibrium bond lengths and binding energies.

Morse [6] approximated the analytical expressions occurring in the Heitler-London treatment with

(6)

where is the molecular binding energy, also called dissociation energy, the equilibrium bond length, and a dimensionless number. Historically, Morse did not introduce the parameter but rather , which is still commonly used today. Moreover, he expressed the energy relative to the ground state, ,

(7)

To ease comparisons between different potentials and to yield , we choose forms similar to Eq. (6)(6) throughout this article.

Morse [6] recognized that the vibrational levels of many diatomic molecules, both homo- and hetero-nuclear, could be described quite accurately when treating , , and as adjustable parameters. The parameter allows the curvature of in the minimum to be adjusted independently from the ratio , however, no additional fine tuning of higher-order derivatives is possible. Given its simplicity, the Morse potential approximates reference data impressively well, as can be appreciated in Figure 3 for Cu₂ and C₂. It even represents the curves obtained for H₂ at different accuracy levels at which the hydrogen molecule is described and substantially better so than cubic fits to the energy minimum, see Figure 4.. Of course, there are also exceptions, where the Morse potential fails, e.g. Cr₂, which is described in Sect. 3.4.

An important result of the Heitler-London treatment is that short-range repulsion between atoms increases (approximately) exponentially with decreasing distance as two atoms approach each other. All highly-accurate potentials use repulsion that does not stray too far away from such a dependence as the pressure of (simple) bulk systems increases approximately exponentially at high compression with pressure for all bonding classes, at least as long as external pressures do not substantially exceed the bulk modulus , which we will get back to in Sect. 6.1.

The simplest generalization of the Morse potential is a double exponential [62] of the form

(8)

which reduces to Morse if the additional parameter is set to . As Abell [62] noted, “this choice is based in part on analytical convenience, but also on the physical grounds that atomic orbitals decay exponentially with “. However, these long-range asymptotics are irrelevant at very small separations between nuclei, as briefly discussed in Sect. 3.2.2.

3.2. Mie/Lennard-Jones and Born-Mayer/Buckingham potentials

The so-called Lennard-Jones potential [8] is the arguably most used potential in molecular simulations. It describes the interactions between entities with closed-electron shells quite well. This includes not only interactions between noble gas atoms, but also the intra- and intermolecular forces between, say, two CH₂ repeat units in polyethylene as long as the two units are not directly connected by a covalent bond. The functional form of the regular LJ potential is

(9)

Here is the binding energy and the equilibrium bond length. The more common way to write the LJ potential is

(10)

where and . However, comparison to other potentials and analytical manipulations are more easily done when starting from Eq. (9)(9) .

The attractive term is the leading-order dispersive interaction. The exponent can be motivated rigorously from perturbation theory, as sketched in Sect. 4.3. In contrast, the exponent from the repulsive term is a mixture of convenience and good luck. It results from the squaring of the term, which was particularly beneficial at times when numerics was not yet done by machines but by humans. This facilitated the work of those who were dismissively called Rechenknechte by theoretical physics professors in Germany in the beginning of the 20ʹth century. The term translates to computing or arithmetic servants and is nowadays used for computers. However, when treating the exponent as a fit parameter in a generalized Lennard-Jones potential

(11)

the exponent generally turns out close to its generic value of .

Replacing the exponent 6 in Eq. (11)(11) with a number yields the Mie potential [80]. The regular LJ potential can therefore be regarded as a limiting case of Mie’s functional form, a so-called 12–6 Mie potential. In contrast, it may not be entirely appropriate to downgrade a Mie potential with arbitrary exponents as a generalized LJ potential, since Mie’s work [80] predated those of Lennard-Jones [8,81] by more than two decades and in fact, it is not clear to us who was the first to use as exponent in the repulsion. Jones, later known as Lennard-Jones, assumed it to be [81] for argon.

Rather than tweaking the exponent of a Mie potential, it is probably more meaningful to replace it with an expression that somehow accounts for the Pauli repulsion between closed electron shells, whose density decreases as an exponential function of the interatomic distance. Slater [82] was the first to identify such a dependence for helium dimers by using a similar approach as Heitler and London [78]. The functional form he identified as being most appropriate for a limited range of interatomic distances was a single exponential, as in the Morse potential. Later, Born and Mayer [83] used the same functional form for the description of repulsion in simple ionic crystals without further theoretical justification but merely by arguing that it fits experimental results better than an inverse power as in Mie repulsion [83]. Buckingham [7] made a similar observation for noble-gas solids. Using exponential repulsion and attraction then leads to the two-body potential

(12)

which is often called Buckingham potential () but also Born-Mayer potential () for two ions of dislike charges, in which case the potential parameters must be constrained to satisfy

(13)

where is the charge of one ion and the vacuum permittivity. Born-Mayer potentials may also contain dispersive interactions in addition to exponential repulsion and Coulomb potentials. For two like ions, there is no bound state and repulsion is dominated by Coulomb interactions at typical interionic distances.

Later, Buckingham [84] found that the repulsion between hydrogen atoms in the lowest triple state and that between helium atoms are better described when replacing the prefactor of the exponential term with an appropriate polynomial. He motivated this with theoretical results obtained using Heitler-London type approaches, which is in line with more recent comparisons [85]. Thus, the term Buckingham potential may also refer to a potential, in which the constant prefactor to the exponential repulsion is replaced with a polynomial in . Van Vleet et al. [86] provide a contemporary discussion on how to construct the polynomials.

3.2.1. Higher-order dispersion and induction

Dispersive interactions are not limited to the leading-order dipole-dipole interactions, yielding the attraction. For example, dipole-quadrupole interactions cause a correction of

(14)

while the quadrupole-quadrupole and the dipole-octopole interactions lead to corrections that asymptotically scale as . Cipcigan et al. [87] recently summarized the hierarchy of dispersive interactions together with those resulting from electrostatic induction.

The most important correction to Coulomb interactions in ionic systems results from the large electrostatic, dipolar polarizability of anions, which was first considered by Rittner [88]. The energy gained by a dipole placed in an external electrostatic field is , while the on-site energy required to create it is , where is the polarizability of the anion. The dipole adjusts itself so that is minimized. This yields a correction of

(15)

to the Born-Mayer potential applied to a heteronuclear, diatomic molecule, where . Assuming bare Coulomb interactions, the polarizability of the cation is assumed to be negligible compared to that of the anion and has the unit of volume. As already recognized by Rittner [88], corrections to the correction in Eq. (15(15) ), arise due to the induction of the cation and the mutual induction of cation and anion. In addition, electrostatic field gradients induce quadrupoles and higher-order derivatives higher-order multipoles. Their incorporation may necessitate damping of the interactions [89,90] or further fine tuning of the on-site interaction between different multipoles beyond linear response to ensure a systematic improvements of the potential.

Rittner [88] found that adding this type of correction allowed him to reproduce important molecular properties of alkali halide molecules while assuming integer charges on the ions. He deemed Pauling’s criterion for the fraction of ionic character of a bond (dipole moment divided by bond length times elementary charge) false, because it fails to include the far-from-negligible polarization deformation on the ions. Madden and Wilson [91] provided much ammunition in favor of Rittner’s conclusion by also considering crystalline structures of many ionic binaries.

It is important to note that the Rittner and related higher-order corrections are not pairwise additive. This is because an additional charge or charge distribution would change the induced dipole and thereby its interaction with the first charge. As a consequence, polarizability in many-atom systems is generally not solved in closed form, as described in more detail in Sect. 4.2.

3.2.2. Damping interactions at short distances

Obviously, higher-order dispersion and induction only matters significantly at small interatomic distances at which point the charge distributions of the interacting atoms start to overlap, which leads to a reduction or damping of the interactions. To prevent artificially large attractions from occurring at small separations, pertinent terms of the potentials are multiplied with damping functions . They are generally constructed such that the overall potential goes linearly to zero with at small , while the damping function quickly approaches unity with increasing . To this end, Tang and Toennies [89] proposed to replace a dispersive interaction scaling with according to

(16)

where is a parameter which depends on the nature of the two interacting atoms and potentially also on the index , while

(17)

being the incomplete Gamma function, which, in the case of a non-negative integer , can be expressed as

(18)

The Tang and Toennies damping function is supposedly the most widely used one.

Not only dispersive but also regular Coulomb interactions between charges or dipoles can be damped, rewind, should be damped with functions like that defined in Eq. (16)(16) to reflect the finite extent of electronic shells and to avoid grossly exaggerated attraction at small separation. In fact, as early as 1962, Dalgarno [92] discussed damping, or, “shielding” as it was called at the time, in the framework of shell potentials accounting for atomic polarizability.

3.2.3. Electrostatic screening at large distances

When atoms approach each other much more closely than their typical nearest-neighbor spacing, as it can happen during high-energy ion bombardment, it is more appropriate to assume the bare Coulomb repulsion, between the nuclei as asymptotic reference rather than the exponential repulsion reflecting the large- asymptotics. Here, are the nuclear charges of the interacting atoms or ions. Due to the propensity of matter to be locally charge neutral, electrostatic interactions are screened at “large” distances, i.e. at distances approaching equilibrium bond lengths, as is also the case, for example, between charged colloids in an electrolyte. To lowest order, the bare Coulomb interaction in a charge-neutral system is screened with an exponential , where is called the screening length. The resulting potential is known as the Yukawa potential [93].

Ziegler, Biersack, and Littmark (ZBL) [94] proposed an empirical screening function with

(19)

where only the screening length depends explicitly on the nuclear charges. ZBL may still be the most used short-range repulsion, despite reported improvements [95], which, however, lack the mentioning of specific functional forms.

In practice, it is necessary to switch between short and long-distance asymptotics to predict repulsion at intermediate . This is commonly done through switching functions as discussed in the appendix of Ref. [96].

3.3. Combining rules for two-body potentials

Bragg [97] observed that the lattice constants of many crystals can be reproduced quite accurately when atomic radii are assigned to individual elements and the assumption is made that two nearest-neighbors touch in ideal crystalline structures. This observation implies no rigorous but an approximate constraint for how (two-body) potentials parameterized for individual elements can be combined to mixed interactions. An arsenal of propositions was made in that regard. The simplest and most wide spread for the Lennard-Jones potential are the Lorentz-Berthelot [98,99] rules. They take the arithmetic mean of the length scale parameter , whereby Lorentz predated Bragg’s observations by 40 years, and, the geometric mean of the binding energies [98,99], expressed here using rather than ,

(20)

(21)

The major flaw of Lorentz-Berthelot is that it generally violates a rigorous result for the mixed dispersion coefficient, which is given further below in Eq. (50)(50) . An apparently reasonable approximation to Eq. (50)(50) for combined dispersion coefficient reads [100]

(22)

where is the polarizability of atom or ion , while the geometric mean merely provides an upper bound for [100]. We abstain from reviewing combining rules not reproducing Eq. (22)(22) or better motivated combining rules, other than Lorentz-Berthelot. Instead, we content ourselves noting that Tang and Toennis [101] found the potential depths and locations of hetero nuclear noble-gas molecules to be well reproduced when using their equations (4)(4) and (5(5) ), which satisfy Eq. (22)(22) from this work, as combining rules. We would expect that reasonable answers can also be obtained for closed-shell (united) atoms when using Eq. (22)(22) in conjunction with the Lorentz rule, from where the prefactor for the repulsion can be deduced.

A better model for repulsion than that used by Mie [80] or Lennard-Jones [8] is the exponential repulsion from Born and Mayer [83], originally Slater [102], which we write here as

(23)

where typically lies in the range of a few to several dozen keV [103]. Values for and pertaining to our selected reference structure can be read off from Table. 1 by associating with and with . Assuming that repulsion originates from he overlap of exponentially decaying charge densities yields the harmonic means for the length scale,

(24)

A standard assumption, similar to the Lorentz-Berthelot rule, is to apply a geometric mean for the energy prefactors, . Abrahamson [103] compiled a list for and for all neutral elements up to the atomic number 105.

An alternative combining rule for the prefactors was proposed by Smith [104] having in mind atoms or ions with noble-gas electronic configuration. He argued that repulsion is an energy excess related to the deformation of the electron density caused by the Fermi principle. He treated that deformation in a way as if each of the two interacting atoms or ions were in contact with an infinitesimally thin but impenetrable wall, which is positioned a distance from atom . The excess energy associated with each side of the wall can then be calculated from the pair repulsion. This nonrigorous, yet plausible model leads to the combining rule

(25)

An identical rule had already been found empirically by Gilbert [105] for alkali halide monomers a few years earlier. Besides proposing his semi-empirical explanation, Smith [104] added to that list various pairs consisting of noble gas atoms, metal ions, and halogen anions, each with a complete shell. Böhm and Ahlrichs [106] found the combining rule for in Eq. (25)(25) to be more accurate than others and further complemented the list of investigated pair repulsions to diatomic molecules beyond those consisting of atoms or ions with closed electron shell. Recent work [86], which is motivated by the analysis of electronic-charge density overlap between two atoms, proposes altered combining rules, in which is the geometric mean of and . In addition, more complicated mixing rules apply than those “derived” by Smith and the prefactors have a certain dependence.

Given its success, it might be in place to further comment on the arguments leading to Eq. (25)(25) . First, it explains quite naturally that repulsion is not necessarily a pair-wise additive quantity. The electron distribution of a first atom is more easily squished by a second atom if there is no third atom already squishing the maltreated electron cloud of the first atom. The poor electrons simply have no more volume of refuge. Thus, the more atoms squeeze against a central atom, the larger the central atom should appear to be. Table 1 reflects this trend for single-component systems not only for the equilibrium lengths but also for . It is particularly strong for the covalently bonded carbon atoms and still quite noticeable for copper. In fact, from own unpublished work on alkali metals, our impression is that the assumption of pairwise additive repulsion is particularly poor for small coordination numbers. Second, potentials adjusting atomic charges on the fly might want to include the effect that this charge has on the atomic size and thus on its repulsive interaction. Neither of these two points appears to have attracted the attention it deserves.

3.4. Funky two-body potentials

Simple, computationally lean potentials generally lead to simple behavior without producing the intricate dynamics of highly viscous liquids or the complex structure of amorphous solids. The phase diagram of Lennard-Jonesium, like that of noble gas atoms other than helium consists of a closed-packed structure at low temperature, a low-viscosity fluid phase in a narrow temperature range at (what would be characteristic for) ambient pressure, and the gas phase at large temperature. While nature manages to produce elements with extended liquid phases at moderate pressure and temperature, like mercury or gallium, they do not become highly viscous. Even binaries often do not form highly viscous fluids and instead either phase separate or crystallize upon cooling from the liquid phase. Notable exceptions to this rule at elevated temperature are SiO₂ and CuZr, which both resist crystallization upon cooling to a significant degree while being describable with relatively simple potentials, e.g. by the silica potential proposed by van Beest, Kramer, and van Santen (BKS) [107] or by the embedded-atom model (EAM) for CuZr [108–110]. Unfortunately, BKS necessitates long-range electrostatic interactions to be evaluated while EAM lacks pairwise additivity, even if it can be computed at the cost of order , like pairwise additive potentials, where is the number of atoms within the cutoff radius and the number of atoms.

To study processes taking place on times exceeding vibrational periods by several decades and to address generic questions typical for disordered solids or highly viscous fluids, it can be advantageous to work with potentials producing the observed macroscopic behavior while being computationally as lean as possible. Such potentials may not be representative of any real material, but, so the hope, produce the correct qualitative behavior for the right reason. In the worst case, they provide a model for what nature could be and thereby allow theories for the fracture of disordered solids [111,112], the supercooling of liquids [113,114] or the rigidity of glasses [115,116] to be tested. Moreover, for many of the frequently qualitative questions to be answered, it can be advisable to sacrifice accuracy in the interactions rather than to make compromises in system size or cooling rate. For example, the anomaly of the specific heat in a bulk-metallic glass former can have serious artifacts when the linear system size is not at least twice the density correlation length [117]. Likewise, determining the proper scaling for the vibrational density of states with frequency in amorphous solids requires the use of large system sizes and long simulation times [118–121] and thereby the use of lean potentials.

The common recipe to construct simple potentials keeping systems from crystallizing quickly is to build frustration into them. Dzugutov [122] achieved this by introducing a hump in the pair potential. Its functional form is given by

(26)

where the parameters , , , , and used for Figure 5 are those from the original work. The most important feature of the Dzugatov potential certainly is that the hump is located near next-nearest neighbors spacings in a closed-packed structure if nearest neighbors settle in the vicinity of the energy minimum. As a consequence, atoms like to adopt local icosahedral structures, which can arrange similar to the order found in quasi crystals lacking long-range order, though large system sizes and small cooling rates may be required to prevent the crystallization into a bcc structure [123].

Figure 5. Two examples of funky pair potentials, , normalized to their binding energy and their bond distance . The solid blue line represents a Dazugutov potential, which was constructed to make a mono-atomic system avoid crystallization. Circles depict experimental results on the pair potential of the chromium dimer, Cr₂. Dashed lines are Morse potential fits to the energy minimum.

For mixtures of particles, frustration can be encoded to some degree through artificial combining rules in Lennard-Jones binaries, as for example, by setting all while choosing , , [113]. This first attempt by Wahnström is quite prone to crystallization as revealed by Kob and Anderson [114]. In an attempt to mimic a Stillinger-Weber potential designed to mimic amorphous NiP [124], they suggested a binary Lennard-Jones potential in which dislike atoms have a large binding energy while being incompatible in size: , , , and , , and [114]. To further reduce the risk of having unnoticed or even worse noticeable crystalline reference phases, continuous distributions of Lennard-Jones radii can be used [115,116]. For a comparison of different glass-forming models, the reader is referred to a recent work by Ninarello, Berthier, and Coslovich [125].

We note in passing that pair potentials with humps similar to the one present in the Dzugutov potential are occasionally used to model mono-atomic bcc phases, in particular that of iron. This is because the third neighbor shell in bcc sits at times the radius of the nearest-neighbor shell and still times that of the second nearest-neighbor shell. Thus the first two shells comfortably fit into the potential well while the third shell has to pay no penalty. In contrast to bcc, the next-nearest neighbor shell in fcc or hcp would be utmost unhappy.

Interestingly, real pair potentials of bcc-forming metals can deviate substantially from the Morse potential, as is the case, for example for the chromium dimer, Cr₂, whose pair potential is included in Figure 5 and which appears to be challenging to compute accurately even with the most advanced post-Hartree Fock methods [126]. Despite the funky shape of “exact” pair potentials of some transition metal dimers, they still do not allow quantitative predictions to be made for elemental bcc metals, because real interactions simply happen not to be pairwise additive. It would lead to a sometimes substantial overestimation of the shear modulus, whereby dislocations are artificially suppressed in crystalline materials (their energy depends predominantly on the shear elastic constants) and thus brittle fracture be enhanced. Pair potentials also lack the important effect that a (metallic) bond between two atoms is strengthened when at least one of the two reduces its coordination as it would happen in a propagating crack. This lack biases material behavior toward brittleness.

Finally, pair potentials with a hump can also be designed to mimic chemical reactions as they occur, for example, during polymeric chemical reactions of step growth [127] and also radical-reaction based chain-growth [128]. It remains a philosophical question if such potentials should be deemed reactive, since the notion of a chemical reaction traditionally requires electron transfer or hybridization changes resulting in altered pair interaction, as is induced, for example between two hydrogen atoms due to the presence of an oxygen atom. Pair potentials can mimic neither one by definition so that they do not classify as reactive, in our opinion. Nonetheless, a well-designed bump potential can reproduce steric effects while reproducing the large energy barrier that needs to be overcome for two monomers to bond, whereby the complex interplay of boundary condition, chain relaxation, and chain growth can be studied.

3.5. Relating potentials and elastic properties

Elastic properties are among the most important properties of solids. Any potential should therefore be tested for its ability to reproduce the elastic tensor of crystalline references. Ideally, they include not only thermodynamically stable but also hypothetical structures or finite stress to enhance transferability. Non-existing reference data can be obtained from DFT or related methods. One of the problems surfacing quickly is that elastic constants are not uniquely defined, except at zero stress and temperature. Different definitions lead to deviations of similar order as the external stress [129–136]. This can matter, for example, under geophysical conditions, so that properly converting between different elastic tensors may be critical, which we will come back to briefly at the end of this section.

To set the stage for many of the calculations throughout this article, we briefly review central aspects of the theory of elasticity, while deriving the Cauchy relations, whose violation has been and remains central to the development of potentials. We attempted our summary to be more condensed than other texts, while being palatable for readers not familiar with the topic.

We assume a deformation of a crystal, in which atomic positions after deformation are . Here, atoms are enumerated by , Cartesian components are represented by Greek letters, and the Einstein summation convention applies. The upper index 0 in denotes an (equilibrium) position in a reference structure, the tensor is the macroscopic displacement gradient, while is a local atomic displacement. The deformation is affine if all vanish. The Eulerian strain tensor is the symmetrized displacement gradient, . From this definition, the squared distance between atoms and , , after an affine deformation can be written as

(27)

where is called the Lagrangian strain tensor [137]. Other strain tensors exist but are not relevant here.

Stress tensor elements are generally defined as the first derivative of a thermodynamic potential w.r.t. a strain tensor element and normalized to the volume of the reference, e.g.

(28)

is called a Cauchy stress. Here, we have normalized energies and volumes per atom (pa). Note that the evaluation of the derivative at only means that the derivative is taken at the reference but not that the stress disappears. Other stress tensor definitions exist, but they only differ when evaluated at nonzero or finite strain. For example, Eq. (28)(28) defines the second Piola-Kirchhoff stress, if evaluated at nonzero .

For systems interacting through (central) pair potentials, can be easily computed. To this end, we first reexpress as so that

(29)

where the prime in indicates the first derivative. For simple lattices under isotropic stress, the summation can be sorted according to neighbor shells, where denotes the nearest-neighbor shell, the next nearest-neighbor shell, and so on. Thus,

(30)

Here is the distance between a “central” atom and an atom in shell and is the second rank shell tensor, whose elements are defined as

(31)

For sufficiently symmetric shell structures, the second-rank shell tensor simply turns out to be , where is the number of atoms in neighbor shell and the spatial dimension [138]. In static equilibrium, the external hydrostatic pressure is nothing but .

Elastic constants are defined as the change of stress with strain, i.e. as second-order derivative of a thermodynamic potential w.r.t. to strain. This implies that results differ at non-zero stress for Eularian and Lagrangian strain tensors, see also Eq. (38)(38) . Moreover, it matters if atomic reference positions, i.e. the Wykhoff positions of atoms in a unit cell are kept fixed or change with strain, as could be expressed by assuming the to depend on strain whenever applicable. Such structural non-affine relaxation can happen when atomic sites lack inversion symmetry, as when shearing a diamond structure, or, as an extreme example, amorphous materials [139].

In-silico, it is an easy matter to constrain all to zero, experimentally less so. Elastic constants defined for are marked with an upper zero. In that case, taking the derivative of the two sides in Eq. (29)(29) simply yields

(32)

By introducing the fourth-rank shell tensor,

(33)

and by using , the elastic tensor reads

(34)

Structural relaxation always reduces the energy so that of diamond would violate the Cauchy relation even more than already evidenced in Figure 2.

The elastic tensor is invariant w.r.t. to any permutation of its indices, including for example, , or, in Voigt notation . It also includes implicitly the recipe for what elastic tensor obeys the Cauchy relations at finite hydrostatic pressure , namely the second-order derivative of the internal energy w.r.t. the Lagrangian strain. However, most simple crystals violate it as is obvious from Figure 2.

Equation (33)(33) allows the elastic tensor of simple crystals obeying pair potentials to be computed in a straightforward fashion, in particular when knowing the fourth-rank shell tensor. Its independent components are compiled in Tab. 2 for selected lattices along with other information useful to quickly compute properties of simple crystals.

Table 2. Useful dimensionless numbers for selected crystal systems using standard orientation of the lattices: simple cubic (sc), face-centered cubic (fcc), body-centered cubic (bcc), diamond cubic (dc), and hexagonal close packed (hcp). Additional elements for hcp are , , and , the negative sign applies to A-layer atoms and the negative for B-layer atoms if B is shifted by w.r.t. A. In diamond, each atom on (0,0,0) and (1/4,1/4,1/4) contribute and , respectively.

sc	6	12		1	1/3	1/3	0
fcc	12	6			1/3	1/6	1/12
bcc	8	6			1/3	1/9	1/9
dc	4	12			1/3	1/9	1/9
hcp	12	6			1/3	5/24	5/72

For short-range potentials only the first two shells contribute substantially to the elastic tensors. Using the information collected in this section so far, we find for fcc, bcc, and sc, using Voigt instead of tensor notation,

(35a)

(35b)

with the effective shell spring constant

(36)

Note that for fcc, bcc and sc, since their primitive unit cell contains only a single atom.

To obtain the relation between different elastic tensor definitions, it is useful to deduce from the definition of the Lagrangian strain that

(37)

Thus, for a second-order derivative of an arbitrary function , it follows that

(38)

The differentiation rule of Eq. (38)(38) also allows one to determine the second-order derivative of the term , where is a constant (external) pressure. Given the relative change of a volume element, and using Voigt notation, can be deduced to be

(39)

with for all , for and both , and else .

Finally, we note that for a solid to be stable, the elastic tensor must be positive definite, which leads to the Born stability criteria [140,141], e.g. , , and in cubic materials. At constant pressure, must be positive definite. For a more general discussion on how imposed boundary conditions affect mechanical stability, we refer to an instructive work by Wang et al. [135], who studied a simple model for gold.

4. Many-body potentials for closed-shell systems

This section focuses on the description of many-body effects in systems composed of atoms and ions in which the constituents can be said to have a closed electron shell. Paradigms are noble gases and alkali metal halides, however, many of the methods and insights conveyed here pertain to a broader context such as valence force fields. The central goal of this section is to highlight the role that atomic dipoles play in the interaction between atoms. These dipoles may originate from quantum mechanical ground-state fluctuations or be induced by an electrostatic field, which arises, for example, on a chlorine atom in rock salt due to a lattice distortion breaking the local inversion symmetry.

Incorporating many-body effects of atomic dipoles—as well as higher-order electrostatic multipoles, which are somewhat meant to be referred to implicitly whenever the term dipole is mentioned—can be achieved in many different ways. For neutral atoms, their effect can be incorporated into potentials in the spirit of the expansion of Eq. (1)(1) . However, as soon as ions linger around, it is more efficient to model the dipoles explicitly. This can be done either by placing a formal dipole on the atom, or, by introducing a Drude model, that is, by coupling the center of mass of the displaced electron shell with a spring to a nucleus, or, to a given interaction site in a molecule. Finally, the dipole can be treated either classically, or, quantum mechanically. The latter appears to be a promising route to model many-body dispersion quite accurately, albeit at a large computational cost. For reasons of completeness, we state that Drude models would better be named after Lorentz, since Drude [142] considered free electrons in a metal, while Lorentz [143] attached (dissipative) springs between electrons and nuclei to model the optical response of bound charges in insulating matter.

4.1. Explicit many-body dispersion

Axilrod and Teller [144] and independently Muto [145] (ATM) extended the second-order perturbative treatment of the quantum dipole fluctuations leading to the leading-order London forces to third-order perturbation theory. The resulting three-body ATM potential reads (for three like atoms)

(40)

where the sum runs over all unique triangles formed by three atoms --. Interior angles of the triangle formed by the three atoms -- are denoted by the Greek letter , where the first index denotes the atom at the apex of the angle. The cosines can be computed straightforwardly from the individual bond lengths, . A frequently used approximation for (mixed) dispersion coefficients is obtained from geometric averages in atomic units, i.e. , where , , are element specific and is 1 Hartree. More rigorous results can be deduced from the Casimir-Ponder integral as described in Refs. [146,147].

The usual course of action when using ATM in conjunction with LJ or Buckingham is to explore improvements on predicted properties. Traditionally, ATM corrections to the binding energy for noble gas crystals like argon are estimated to be of order 10% [49]. For very polarizable systems, they can be much larger. Von Lilienfeld and Tkatchenko [146] found them to account for 50% of the binding between two graphene sheets. In addition, ATM interactions can correct for systematic deficiencies in the elastic and vibrational properties that cannot be overcome in the realm of pair potentials [148]. We see it as beneficial to analyze as an independent quantity in terms of a lattice ATM constant, which can be defined in analogy to the Madelung constant as

(41)

We obtain for an ideal fcc crystal from which one sixth can be assigned to each atom. For hcp, we obtain , which can topple the balance in favor of fcc. The differences between these two numbers is small, yet, larger than between the results for , i.e. , for fcc, versus, , for hcp. In order for to be larger for hcp than for fcc, interactions beyond next-nearest neighbors must be included. Unfortunately, we did not manage to find recent literature results on the simple-to-compute ATM lattice sums beyond relatively rough, initial estimates by Axilrod [149]. Nonetheless, different authors [150,151] have previously concluded that three-body dispersion favors fcc over hcp, although the corrections due to nuclear quantum fluctuations appears to be substantially more important in that regard [152].

4.2. Classical polarizable potentials

The electrostatic field on a central atom or molecule is produced by other charges, dipoles, or higher-order multipoles. This leads to a hierarchy of inductive interactions [90], which we have already touched upon in Sect. 3.2.1. It contains the attraction between a charge and an induced dipole. Next in line is the (asymptotic) attraction between a permanent dipole and an induced dipole or that between a charge and an induced quadrupole. Unfortunately, the interactions related to induced dipoles or induced higher-order multipoles, which will be ignored for the most part in the following, are not pairwise additive. Two like charges placed at from an atom will not induce a dipole on that atom and thus not lead to twice the energy gained if only one of the two charges were present.

Since adjacent existing (static, molecular) dipoles try to align themselves to an external electrostatic field in an attempt to minimize the potential energy, induced dipoles tend to be parallel to pre-existing static dipoles. This effect is well known to increase the mean dipole moment of water molecules from 1.85 D in the gas phase to approximately 2.7 D in the liquid, whereby dipole-dipole interactions essentially double.

The polarizability of a homogeneous medium having cubic or higher symmetry is stated in terms of its dielectric constant , which implicitly reflects the feedback that dipoles have on each other. In atomic or molecular systems, can be well estimated through the Clausius-Mossotti relation [153]

(42)

where is the number density of species and its (orientationally averaged) polarizability. Once the right-hand side of Eq. (42)(42) is greater or equal one, is no longer a finite positive number so that the model has reached its physically meaningful limit. The dipoles “want” to grow ad infinitum, at which point the system becomes metallic. The underlying polarizable potential needs to be extended to either suppress this so-called polarization catastrophe, e.g. by shielding the Coulomb interaction at small distances [89], or, to actually allow the system to become conducting. In advanced shell potentials being a compromise between polarizable and charge-transfer models, this can be achieved by not constraining an electron cloud to one particular atom [154].

When modeling charged or polar systems by atomistic means, one would certainly want to reproduce the dielectric response function of a medium correctly, e.g. when simulating the condensation of water on a surface [155], the Helmholtz double layer on an electrode [156,157], or the damped Coulomb interaction of charged colloids in water. For the simulation of simple geometries and homogeneous media, it may be possible to achieve this with effective potentials [158] or with the concept of induced mirror charges [159]. However, many important problems lack the symmetry or isotropy required to pursue such approaches so that the dielectric response has to be solved on the fly. This can be achieved with molecular approaches encoding the polarizability into the potential energy surface.

One possibility is to induce ideal dipoles or higher-order multipoles on atoms or specific interaction sites in a molecule [160,161]. This is done by adding an energy contribution for each inducible multipole . The usual interatomic potential is the one that minimizes , e.g.

(43)

where

(44)

with respect to the dipoles. In Eq. (44(44) ), contributions to the short-range (sr) interaction are separated from those due to Coulomb (C) interactions, which may contain contributions from static dipoles (for molecular rather than atomic simulation) in addition to those from point charges. Since, Eq. (44)(44) is a second-order polynomial in the induced dipoles, a well-defined minimum exists as long as the Hessian related to the dipoles is positive definite. A brute-force inversion of the Hessian to yield the exact minimum is generally inadvisable. Alternatives are conventional minimization techniques, such as those based on conjugate gradients [162] or extended Lagrangians [163,164]. In the latter case, the dipoles, or other adjustable degrees of freedom such as those describing the shape of the periodically repeated system [163] or prefactors to electronic wavefunctions [164], are assigned an inertia and propagated and relaxed along with the atomic coordinates. It is beyond the scope of this review to discuss the pros and cons of extended Lagrangians in detail. It suffices to say that their implementation is relatively simple. However, they introduce an effective delay on the molecular dynamics [165]. In addition, relaxation to the energy minimum can be slow when the Hessian has strongly differing eigenvalues, which automatically happens for large . Then, the dipolar response functions are “stiff” at small length scales but “soft” at the continuum scale.

One disadvantage of placing ideal dipoles on atoms or ions is that the evaluation of Coulomb interactions is substantially complicated, even if solutions exist to include their effect into the (fast) Ewald summation [166]. Another deficiency is that higher-order multipoles are ignored that are generated when an electron cloud displaces with respect to a nucleus. These drawbacks can be remedied with Drude particles, or “Drudes”, in which a fixed Drude charge is coupled harmonically with a spring constant to an atom or ion, to which a charge is added [167–169]. The last summand on the r.h.s. of Eq. (44)(44) must then be replaced with , where is the displacement of one Drude charge w.r.t. the atom or interaction site that it is bonded to. On-site Coulomb interactions between all charges on the Drude, which despite similar spelling is not to be confused with The Dude from The Big Lebowski, must be switched off. Minimization of the total energy w.r.t. to the Drude displacements can be done in a similar fashion as for ideal dipoles.

In order for the Drude to reproduce the correct polarizability, the relation

(45)

must be obeyed. The limiting case of ideal dipoles is obtained for while keeping constant. With an appropriate choice of and sign for , the correct quadrupolar induction of a spherically symmetric Drude can be matched to a homogeneous field. Traditionally, the three independent parameters of Drude oscillators, which in addition to and are also assigned an inertia or mass , were chosen to best match the frequency dependence of the dielectric constant at high frequencies [167]. Ideally, those values would be close to results obtained from a fit to accurate reference data on forces and energies, for which and possibly even are treated as adjustable parameters. Large discrepancies between different parametrization procedures should probably be seen as a sign that something is missing or wrong with the potential.

A stringent test for the correctness but also for the relevance of dipolar polarizability can be obtained when comparing force-field-based infrared absorption spectra to reliable reference data. For example, in the case of amorphous silica, obtaining the correct number, positions, and intensities of peaks requires the polarizability of the oxygen atoms to be included [170] (Figure 6a). Accounting for electrostatic induction also appears necessary to accurately reproduce bond-angle distributions in silica and related systems [171]. While global bond-angle histograms may seem fine for pair potentials, deviations between symmetry-specific angles in crystalline structures from the ideal tetrahedral bond angle have the wrong sign for the most commonly used SiO₂ pair potential [107] but the correct sign and magnitude for the polarizable Tangney-Scandolo potential [171,172].

Figure 6. Influence of polarizability on properties of SiO₂. (a) Infrared absorption spectrum of amorphous SiO₂ (silica) computed one time with a rigid, fixed-charge ion model (dashed line) and one time with a dipol-polarizable anion (solid line). The polarizable model matches qualitatively the experimental results shown in the inset. (b) ratio of crystalline SiO₂ (quartz) computed with the rigid ion potentials of van Beest, Kramer and van Santen (BKS) [107], the fluctuating-charge potential of Demiralp, Çaǧin and Goddard after modifying parameters (mDCG) [173], and the polarizable force field of Tangney and Scandolo (TS) [171]. Only the polarizable force field reproduces the anomaly at the --quartz transition with temperature . (a) is reprinted with permission from Wilson, M., Madden, P. A., Hemmati, M., and Angell, C. A. Phys. Rev. Lett. 77, 4023 (1996) (Ref. [170]). Copyright (1996) by the American Physical Society. (b) is reprinted from Herzbach, D., Binder, K. and M. H. Müser, J. Chem. Phys. 123, 124711 (2005) (Ref. [172]). with the permission of AIP Publishing.

The incorporation of dipolar polarizability is particularly crucial for structures in which highly polarizable (united) atoms or ions, most notably anions, are located in sites deviating strongly from inversion symmetry. This concerns in particular oxygen whenever it is two-coordinated, as it is in water, but also in low-temperature tetrahedral network formers like silica. If competing local orders exist, neglecting the polarizability of anions will artificially favor the order with the smaller deviation from inversion symmetry, because the anion is shorn of its ability to reduce its energy in the larger electrostatic field of the symmetry-broken phase.

Including dipolar polarizability can also be required to describe macroscopic structural changes occurring during the transition between high- and low-symmetry phases in a qualitatively correct fashion, as is the case for the transition in quartz. Using pair potentials, the ratio shows no anomaly as a function of temperature at the transition temperature. For charge-transfer potentials, the slope of the ratio changes at the transition temperature, while it is discontinuous, as in experiment, when dipolar polarizability on oxygen is added to the pair potential [172] (Figure 6b). This discussion may show that the analysis of phase transformations can benefit the validation of potentials, most notably polarizable potentials, even if applications will evolve primarily around liquid or amorphous phases.

An affine deformation of a highly symmetric ionic crystal, as for example, rocksalt (B1) or cesium chloride (B2), does not lead to an electrostatic field on a central lattice position created by the atoms residing on other positions. Thus, the central atom does not develop a dipole. However, a non-isotropic deformation, e.g. a uniaxial strain, can reduce the cubic symmetry, whereby the charge density on the anion may develop a quadrupole moment. Consequently, unlike dipolar polarizability, quadrupolar polarizability can induce a Cauchy violation in the said crystals, which is why it can be said to be more important than dipolar polarizability, at least from a continuum-mechanics point of view. Madden, Wilson, and coworkers demonstrated that inclusion of quadrupolar polarizability can substantially increase the agreement between experimental and simulation results, at least in the case of simple ionic systems like AgCl [161] or MgO [174]. An important aspect of their work [174] is that they managed to gauge adjustable parameters independently from one another, even those describing a charge-density change on an anion due to electrostatic polarization from those being caused by short-range repulsion.

4.3. Quantum-Drude oscillators

The starting point for the description of dispersive interactions are two, or more, isolated atoms, described by their free-atom Hamiltonians . The atoms interact through their dipoles via plus potentially through higher-order multipoles. Second-order perturbation theory then leads to a pairwise-additive interactions between the atoms.

Asymptotic interactions are also obtained when replacing the atomic Hamiltonians with quantized Drude oscillators [90,175,176], as will be shown here below. Their three, free parameters per atom, , , and can be chosen to match the correct polarizability and leading-order dispersive coefficient , which allows dispersive and inductive interactions to be treated on a common footing. The justified hope is that higher-order and many-body dispersive interactions can be mimicked reasonably well simultaneously by properly selecting the third remaining parameter [87]. If the nucleus is also quantized, must be replaced with a reduced mass .

The Hamiltonians of two identical, three-dimensional quantum Drudes can be decoupled into three pairs with

(46)

The parameters (, no summation convention) take the values or depending on whether the displacements are orthogonal or parallel to . The oscillators can be decoupled further through the transformation leading to three isolated oscillators pairs described by Hamiltonians of the form

(47)

In leading order, their combined excess ground-state energy w.r.t. that of two uncoupled oscillators, with per free quantum Drude pair is

(48)

Adding up the three Drude pairs, using , and repeating the entire exercise for potentially dislike Drudes, and , then yields,

(49)

which is the same combination rule as Eq. (22)(22) . Using the true atomic reference Hamiltonians and placing the dimer on the -axis, the correct mixed dispersion coefficient reads [177]

(50)

where is the matrix element of an atom and where the primed sum indicates that at least one of the two quantum numbers and must differ from zero.

Typical parametrizations of quantum Drudes [90] yield of order 0.5 H to 1 H for hydrogen and noble gas atoms as well as for small, closed-shell molecules, in between 0.7 and 1.4, and from 0.1 to 0.6. Since a regular quantum Drude is fully defined by three parameters, many terms related to higher-order polarizability (, dipole, 2 = quadrupole, 3 = octupole) are constrained to take fixed ratios. Some of them are given by [87]

(51a)

(51b)

(51c)

Figure 7 reveals that these ratios tend to reproduce experimental results within 10% error.

Figure 7. Experimental values for the left-hand sides of Eq. (51). Reprinted from J. Comput. Phys., 326, Cipcigan, F. S., Sokhan V. P., Crain, J., and Martyna, G. J., Electronic coarse graining enhances the predictive power of molecular simulation allowing challenges in water physics to be addressed, 222–233, Copyright (2016), with permission from Elsevier.

Once a quantum Drude model is defined, various strategies exist to deduce the energy and interatomic forces following from it. Brute-force diagonalization is not advisable for reasons of computational efficiency but also because the model is no longer quadratic in the Drude displacements, as soon as they seize to be negligible compared to interatomic distances, i.e. at typical condensed-phase, nearest-neighbor spacings. Moreover, any solution strategy should allow for the possibility to replace the harmonic coupling with one in the spirit of finitely-extensible non-linear coupling, as to potentially match simultaneously more than just one inductive or dispersive coefficient in addition to and .

A common strategy to handle quantum Drudes is the use of path-integral techniques [178–181]. The trouble with quantum Drudes is that the number of system replicas to be simulated is of order one Rydberg divided by thermal energy, which is a few hundred at room temperature. One can get away with slightly fewer replicas with diffusion Monte Carlo [182], which is closely related to path integrals, in conjunction with a diagrammatic expansion of interacting Drudes [183]. Unfortunately, improvements do not seem such that numerical costs fall below an approach, in which induction is handled classically and ATM interactions are included explicitly. This turns the quantum Drude oscillator approach into something like a sleeping beauty, which could be awaken if a handsome prince managed to cut down the computational overhead w.r.t. classical polarizable potentials plus two-body like short-range potentials to a factor of order ten or less.

4.4. Shell models and many-body repulsion

The interaction of atoms results in the deformation of (closed) electronic shells beyond the induction of electric multipoles through electrostatic fields, reflected, for example, by an approximately rigid translation of the valence shell with respect to the remaining ion. A reduction of the Coulomb interaction between charges and dipoles due to the overlap of electronic shells can be described with damping functions, which were introduced in Sect. 3.2.2. The electrostatic field, its gradients but also the Fermi principle and hence repulsion deforms the electronic density compared to the superposition of the free-ion or free-atom references. These deformations were first approximated as being spherical in the so-called breathing-shell model [184] and described on a common footing with dipolar polarizability [167]. Later, non-spherical shell deformations were also considered with deformable-shell [185], deformable-ion [161], or distortable-ion [174] potentials. Variables describing the state of deformation are the ionic radius, or rather, the deviation from the radius of the free ion, plus variables describing the shape of the deformation. In atomistic simulations, the shell variables, including the induced dipoles, are assumed to minimize the potential energy.

The starting point of much of the original literature [48,184,185] assumes a harmonic expansion of the potential energy in terms of small variables, i.e. displacements, dipoles, and breathing modes. However, for general situations, the potentials are better cast in terms of implicit many-body potentials, the ingredients of which are the repulsive interactions and on-site energy penalties for the deformation of the shell with respect to free ions. In such models, the short-range repulsion can be expressed, for example, in a breathing-shell potential, through an expression of the type [174,186,187]

(52)

Here, , , , and are constant coefficients, while the minimize . The first summand on the r.h.s. of Eq. (52)(52) reflects the on-site coupling in a heuristic fashion, i.e. it can be parameterized to yield the correct shell stiffness at small while suppressing embarrassing ion shrinkage at large compression.

Accounting for anisotropic shell deformation requires tensors of rank one and higher to be included as arguments in the functions appearing on the r.h.s. of Eq. (52)(52) such that they reduce to a scalar in a meaningful way, i.e. without violating the isotropy of space. Anisotropic shell distortions seem to be required to simultaneously improve phonon spectra and the binding energy differences of competing ionic structures compared to rigid-shell potentials merely reflecting dipolar polarizability [184,188]. We note in passing that aspherical atoms modeled with potentials assuming fixed bonding topography has been pioneered by Price, Stone, and coworkers [189]. Deviations from spherical symmetry then stems predominantly from intramolecular interactions so that shape parameters can be treated as fixed relative to a molecular coordinate system.

A central motivation for introducing flexible shell models was to reproduce the experimentally observed violation of the Cauchy relation in simple salts and their phonon dispersion at high symmetry points [48,184]. Here, it might be worth observing that the Cauchy pressures tend to be minor for alkali halides but quite substantial for alkaline earth oxides like MgO. One central difference between these two classes of simple salts is that a free, singly-charged halogen anion is stable while a free, doubly charged chalcogen anion is not. Thus, assuming doubly charged ionic references without accounting for charge transfer upon a change in bond length could be difficult to justify.

Interestingly, charge-transfer and breathing-shell models have a similar effect on the Cauchy pressure. When the application of a stress in -direction induces ions to shrink in all spatial directions, the external pressure required to maintain the crystal shape in the or direction is reduced, which implies that is reduced compared to a rigid-shell model with fixed or zero . In contrast, when shearing, for example, a rock-salt structure, the nearest-neighbor bond lengths only change in order so that the shear modulus is the same as that of a rigid-shell potential. In other words, , is also unavoidably obtained for charge-transfer potentials applied to rocksalt crystals, irrespective of whether neutral atoms or ions are taken as reference [190]. Figure 8 pictures the just-made arguments for shell models schematically and explains qualitatively how anisotropic shell models can lead to a Cauchy pressure of either sign. Closed form expressions for the elastic constants in breathing-shell models are summarized in Ref. [191], which also contains a critical comparison of different breathing-shell potentials.

Figure 8. Deformation of an electron shell in response to a deformation of the macroscopic body. The electron shell is schematically depicted by the central circle. (a) During shear deformation, the displacement does not lead to a deformation of the shell. There is hence no influence of shell deformation on the shear modulus . (b) If the shell deforms isotropically, then an inwards motion of the black atoms will lead to an inwards force on the white atoms, implying . (c) For anisotropic shell deformation, an inwards motion of the black atoms will lead to an outwards force on the white atoms, leading to . Adapted from , Sangster, M., Interionic potentials and force constant models for rocksalt structure crystals, 355–363, Copyright (1973), with permission from Elsevier (Ref. [192]).

Ion shrinkage under compressive stress without charge transfer and the unavoidable sublinear scaling of repulsion with the coordination number in breathing-shell models strikes us as implausible. It certainly violates the superlinear scaling of repulsion with obtained in a first-order DFT-based perturbation theory, in which the kinetic energy density of the electrons increases with , where is the electronic density. Despite this argument being merely qualitative, it should reflect the proper trend as the trend is very clear. Our own quick and dirty analysis of the EOS of MgO in the NaCl and CsCl structures using exponential repulsion plus Coulomb interaction and the respective Madelung constants makes us believe that trends are complicated. Although the compressive part of the EOS can be described quite well with the most simple Born Mayer potential, the fitted charges turn out close to unity in both cases while Bader analysis finds more intuitive charges close to . When using the parameters deduced for the NaCl and CsCl-structure for the dimer, differences in the Coulomb energy turn out larger than for the repulsion. At the same time, the Bader charges barely change with the lattice constants in contradiction to our suspicion that breathing-shell models reflect charge-transfer effects to low order. Nonetheless, short-range repulsion certainly induces non-spherical distortions leading to electrostatic multipoles, which then must be included in the overall electrostatic interactions between atoms [161].

5. Many-body potentials for open-shell systems

By definition, open-shell atoms have an incomplete valence shell. Bonding between them occurs through the formation of either covalent, or, in condensed phases, also through metallic bonds. Interactions between open-shell atoms cannot be faithfully described with pair potentials without producing false trends, some of which are discussed in Sect. 2.3. The arguably most important many-body effect in open-shell systems is the weakening of a bond due to the presence of additional atoms. An extreme example is the onset of repulsion between hydrogen atoms, which attract each other as free radicals but repel each other in the presence of oxygen. A more subtle effect is the contraction of layers near unpassivated metal surfaces, where missing neighbors of the atoms in the surface layer strengthen their interaction with atoms in the layer underneath.

Figure 9 reveals this weakening for copper and carbon in different crystalline structures, including hypothetical structures. The energy per atom as a function of nearest-neighbor bond length (Figure 9a – copper, Figure 9e – carbon) was obtained from the materials scientist’s favorite electronic structure method, DFT [193,194] within the local density approximation [194] using projector-augmented waves [195]. The minimum of these curves is the cohesive energy of the crystal , where is the equilibrium bond length shown in Figure 9b for copper and in Figure 9f for carbon. is approximately a logarithmic function of the coordination number for both copper and carbon. The energy per bond in the equilibrium configuration, , clearly decreases with and much more so for carbon (once ) than for copper. For both elements, the decrease of is approximately algebraic in , albeit with a steeper power law for carbon than for copper.

Figure 9. (a,e) Energy per atom in selected crystal structures as a function of the nearest-neighbor spacing . Solid lines are fits of the data to a standard Morse potential assuming nearest-neighbor interactions. (b,f) Equilibrium bond length as a function of coordination number . (c,g) Equilibrium binding energy per bond as a function of . (d,h) Equilibrium bond length as a function of equilibrium binding energy per bond (Pauling plot). Crystal structures range from 1D to 3D include hypothetical ones and are abbreviated as follows: dimer/molecule (mol), chain (ch), graphene (gra), triangular lattice (tri), kagome lattice (ka), square lattice (sq), diamond cubic (dc), simple cubic (sc), body-centered cubic (bcc) and face-centered cubic (fcc). The top row shows results for copper, while the bottom row shows results for carbon. and are the binding energy and bond length of the dimer. Data, which were produced using consistent DFT-LDA methodology, are merely meant to convey trends.

The results presented in Figure 9 imply that pair potentials cannot describe the energetics of these systems, although each individual at fixed can be reproduced quite accurately using a regular Morse potential as evidenced by the solid lines in panels (a) and (e) of Figure 9. As Abell [62] emphasized, environment-dependent parameters are needed. Various families of many-body potentials for open-shell systems have been developed over the years to reflect these trends and, of course, more subtleties. In this section, we will outline their functional forms and the properties of the solids that they describe. We will also motivate some of them from quantum-mechanical considerations.

A commonality of the most popular, simple open-shell potentials, namely EAM, second-moment tight-binding expansion (TB2M), Stillinger-Weber (SW), and Tersoff, which will all be introduced in this section, is that their total potential energy can be written as a cluster functional [19]

(53)

where is a pairwise additive interaction—to be distinguished from the pair potential—and is an effective pair interaction, which, however, depends on a bond-environment variable , which is not necessarily symmetric in the indices and may contain a non-zero on-site term . Its specific interpretation changes from one potential class to the next but it always depends on a sum over third atoms through

(54)

In a limiting case, simply corresponds to the third term, , of the expansion given in Eq. (2)(2) . However, can also implicitly include higher-order terms. The way how to express the various potentials in the generic form of Eq. (53)(53) and (54(54) ) are detailed in Sect. 5.6. Those summaries may be useful to design prototype templates for functions or to work out elastic properties falling into this general category.

5.1. Density-dependent potentials

The alkali metals can be seen as a polar opposite to closed-shell systems, because their electrons are delocalized to the extent that their dispersion relation is close to that of free electrons up to the Fermi wave vector , albeit with an effective mass. This is because there is only one electron per atom in the conduction band so that only its minimum is sampled. For this reason, early attempts of describing bonding in (alkali) metals assume the Jellium model [75,196,197], which expresses the energy of an electron gas as a function of its density , where in the case of (neutral) alkali metals the number of atoms is identical to that of valence electrons. The Jellium model is the DFT-based approach to matter, in which the main effect of ions is assumed to provide a charge-neutralizing background to the electrons without explicitly accounting for the relative positions between ions. Interestingly, the Jellium model leads to densities and bulk moduli that are relatively close to those of crystals formed by light alkali metals, even when using the true electronic mass [196]. However, Jellium does not stiffen sufficiently much with increasing density when the Coulomb repulsion between the ionic cores starts to matter.

The potential of the form [198,199],

(55)

can certainly be motivated from the Jellium model. In this potential, which could be coined a global-glue potential, the Jellium model is augmented with a pair-wise additive repulsion, e.g. in the form of damped Coulomb interactions. In principle, not only the pair repulsion but also can be made a function of the element, or, elements under consideration. While Eq. (55)(55) can be generalized to non-homogeneous phases and alloys, e.g. via the construction of Wigner-Seitz cells and the request of local charge neutrality, surfaces would cause problems, since surface atoms “own” excessive volume, which makes a local density difficult to define. Thus, it is clear that potentials footing on Eq. (55(55) ), or, simple generalizations thereof, are not practicable. It is yet utmost instructive, to study their properties and further reasons for their failure.

Assuming interactions past the nearest-neighbor shell to be screened, see Sect. 9, severe restrictions for the elastic constants of simple solids ensue. The glue effectively acts as a positive, external pressure, , squeezing the ions together. Combining Eqs. (35) and (39(39) ) then allows the elastic tensor elements of elementary fcc and bcc crystals at zero external stress to be approximated (assuming a smoothly changing with density) with

(56)

where would be an element-specific constant. Equation (56)(56) predicts . Check mark for simple metals! For fcc, the sum rule can be read off. It is obeyed reasonably well by our favorite metal, and other fcc metals, , , , , but less so by the favorite metals of our wives and .

Unfortunately, bcc turns out mechanically unstable, because the pertinent elastic tensor is not positive definite owing to for positive glue pressures. This inequality can be easily understood. The global glue attempts to increase density. This is done most effectively when nuclei keep a maximum mutual distance, which the most closed-packed structures achieve the best. But why then do alkali metals condense into bcc? Before addressing that question, we note that tends to be quite close to unity for simple bcc alkaline metals but not for bcc metals with an approximately half-filled -shell. Specifically, within 2% for the alkali metals from lithium to cesium, while transition metals assume much larger ratios. For example, 1.95 (V), 4.37 (Cr), 2.60 (W), and 1.76 (Fe). Thus, density-dependent potentials fail distinctly more for transition metals than for alkali metals. Also, most alkali metals transform into fcc at relative moderate pressures less than 10 GPa, thus, check mark again, before they actually undergo a series of additional phase transformations at even higher pressures [200]. Thus, despite being highly flawed, Eq. (55)(55) contains a few elements of truth for alkali metals. However, the existence of non-closed-packed equilibrium structures at large compression is also a clear indication that the assumption of pairwise repulsion may not be particularly accurate.

In a seminal work, in which the term electron correlation appears to have been introduced, Wigner [196] demonstrated that free electrons at small density can lower their energy compared to the cheapest, simple constant-density solution of fermions, which is spanned by the product of two Slater determinants, one for spin up and one for spin down, in which all -states are filled up to . At electron densities characteristic for alkali metals and using second-order perturbation theory, he found them to condense into a bcc crystal. One possible conclusion is that a reliable potential for alkali metals should implicitly reflect higher-order gradients in the electronic charge density. Stabilizing the right phase for the wrong reason, e.g. with potentials merely depending on (estimates for) the local electronic charge density and/or by tweaking cut-off distances or functional dependencies, as done, for example by some funky potentials in Sect. 3.4, cannot lead to accurate, transferable potentials.

5.2. Embedded-atom method and second-moment tight-binding potentials

A simple generalization of Eq. (55)(55) is

(57)

where would be an estimate for the electronic charge density near atom . Equation (57)(57) forms the basis of a class of potentials, which has two names, embedded-atom method (EAM) and tight-binding (TB) approximated to second-moment (2M) expansion. It has been historically pursued by two communities. The TB2M and EAM camps differ predominantly in how they motivate and later generalize for alloys or encode directional bonding, which will be sketched further below. Yet, neither camp keeps the interpretation of as the charge density of all valence or conduction-band electrons at or averaged over the vicinity of and assumes instead

(58)

which excludes the contribution of the valence shell of the atom itself. The function is the square of a bond-integral in the TB interpretation and the charge density from atom seen by atom for the EAM camp. However, as Finnis and Sinclair [47] rightfully note in their pioneering study: “the consequences of the form of the model […] does not depend on the physical interpretation.”

Ducastelle [46] appears to have been first to suggested Eq. (57)(57) building on earlier work by Cyrot-Lackmann [201] and Friedel [202] on TB2M approaches to metals. He also proposed what could be seen as the most generic functional form of an EAM or TB2M potential, simple exponential repulsion and furthermore

(59a)

(59b)

For alloys, and depend on the atom type. Moreover, atom (EAM) or bond (TB2M) specific prefactors must be added to each summand on the r.h.s. of Eq. (59b)(59b) . The central arguments leading to the square-root dependence of are sketched in Sect. 5.8. The main reason for the exponential dependence is, as always, the exponential (Slater) type of the atomic orbitals, even though the full TB bond integrals have polynomial prefactors, see, e.g. Ref. [203]. This functional form was optimized for a wide variety of metals and alloys by Cleri and Rosato [204].

Recast into our generalized Morse form, Eq. (8(8) ), the potential energy expression becomes

(60)

We would like to argue that the usefulness of this expression for rationalizing trends in metallic bonding (see Sect. 6) is on par with that of the Lennard-Jones potential for noble gas atoms, which is why we feel that calling systems described by this potential Ducastellium is as appropriate as using the well-established term Lennard-Jonesium.

Historically, the EAM/TB2M potential described by Eq. (57)(57) was independently discovered several times after Ducastelle’s work. For example, Gupta [64] showed that it reproduces the lattice relaxation of metals near surfaces. Tománek, Mukherjee, and Bennemann revealed its appropriateness to describe the energetics of small metal clusters [205] as well as surface and vacancy energies of transition metals [206]. In thermodynamics, a formulation identical to EAM/TB2M has been employed by Pagonabarraga and Frenkel for coarse-grained particle dynamics calculations, where becomes a free-energy that is adjusted to the equation of state of a liquid [207,208].

Equation (57)(57) can also be motivated from the quasi atom theory of Stott and Zaremba [209] and from the effective-medium approximation theory of Nørskov and Lang [210,211]. In their approaches, the energy gained when embedding an atom into a given site , e.g. an interstitial or vacancy site, is argued to be a functional of the electronic density at the embedding site and thereby to depend in leading order on the charge density that exists at before the atom is embedded. Daw and Baskes [9,45] built on these ideas and approximated the charge density at the embedding site as a superposition of the charge density from neighboring atoms. They coined the term embedded-atom method and provided important insight as to what extent many-body terms, as well as hydrogen, affect materials behavior including ductility. Ercolessi, Parrinello and Tosatti renamed EAM into glue potential [212–214] while claiming that their specific realization of an EAM potential “accounts for all known lattice properties of Au” [212].

Modern EAM potentials use more complex functional forms than the simple exponentials of Ducastelle. Those are out of scope for the discussion in this review and the interested reader is referred to the original literature, see for example works by Mishin and coworkers [215–219]. In the spirit of early works by Foiles, Daw, and Baskes [220], parameterizations of such EAM potentials are distributed using tabulated data for , and the pair-potential .¹ Care has to be taken when interpolating between these data points in the final implementation of a potential, as the choice of spline order is not to be cast aside as a technicality, because it affects the properties of the potential [221].

In our own experience on copper [222], the simple form first used by Ducastelle [46] and Gupta [64] performs best when testing structures with coordination numbers varying systematically from to . The reason for this trend might be, as so often, that simple, physically well motivated functional forms are less prone to overfitting than elaborate functions having been tweaked to enforce right numbers for selected properties in one or few structures. To this we wish to add that we highly doubt the pair-additive repulsion in open-shell systems to be an accurate approximation. Thus, there is little reason to fiddle around with the 4ʹth digit of an embedding function, if “bond-order effects” on repulsive forces lead to errors in the second digit.

5.3. Modified-embedded-atom-method potentials

A deficiency of EAM potentials is their generic preference for closed-packed structures. While funky parametrizations allow bcc to be stabilized, we refer yet again to Wigner’s work on alkali metals in Sect. 5.1 and repeat our claim that stabilizing the right phase for the wrong reason cannot yield a robust potential. To better encode directional bonding in EAM, Baskes [223,224], having silicon in mind, suggested to augment the computation of the electronic density with an angular three-body term. Baskes already hinted [224] that his explicit expression can be interpreted as a dependence of the embedding function on the gradient and higher-order spatial derivatives of the density, i.e. to replace with , where is short-hand notation for . Shortly after, Baskes [225,226] suggested to replace (the estimate for) the embedding density according to with

(61)

the first four satisfying

(62a)

(62b)

(62c)

(62d)

where are weighting coefficients, is a unit vector parallel to bond -, and Einstein summation convention is implied on the Cartesian indices for the squared quantities. Here, the reflect generalized, partial background electron density from atom seen by atom . Through this generalization, semi-explicit angular dependencies result, as for example, through the squaring of the term. It leads to a summand of the form , which is proportional to . With each higher-order term, higher-order sinusoidal dependencies on bonding angles result. An appealing aspect of 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

(63)

where the upper index on the partial embedding charge densities was omitted. clearly has a minimum at the preferred bond angle , since is the cosine of the angle --. 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 with so that after executing the restricted double sum correctly

(64)

Assigning an energy penalty on 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 and correct bond angles can be favored over closed-packed structures. It is clear that Eq. (65)(65) 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 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) [222]. For example, if , then in addition to there can be a contribution of , which, in principle, can be linearly independent of the former, rather than to be constrained to , which can be loosely associated with the term of Eq. (62c)(62c) . Yet another invariant, which can be constructed with four Cartesian indices and a square in the densities, would be . Thus, to flesh out our second point of criticism, the and terms in Eq. (62c)(62c) are far from complete and it is not clear if the most relevant ones have been considered.

The “expansion” of Eq. (62c)(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 can depend, in principle, on . 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 beyond screening, e.g. an environment-dependence of . 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 [227] 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

(65)

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

(66)

with for the tetrahedral angle. Besides this angle, the SW potential has eight independently adjustable parameters: , , , , , for the pair potential, and in addition , 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 [227,229,230]. The latter goes hand in hand with a coordination change from in the diamond structure to in the liquid [231].

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 [232]. 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 , which explains why SW is occasionally used for the united-atom modeling of water [233]. Using small , SW can describe condensed phases of tin, including that condensing in the -tin structure [234], which can be seen as a compromise between tetrahedral and close packing.

Of course, as Stillinger and Weber [227] admitted, their model has quantitative deficiencies. They were revealed most clearly by Biswas and Hamann [235], 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. [236,237]. 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 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 . Such environment-dependent interaction potentials (EDIPs) were first developed for silicon [238,239] and later for other elements including carbon [240,241]. 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, sp²-hybridized) and 4-fold (diamond-like, sp³-hydridized) atoms in amorphous carbon obtained from liquid quenches at varying densities [241].

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

(67)

where are expansion coefficients, Legendre polynomials, and the functions, which can be defined to reproduce the exact or any general three-body potential. When assuming them to factorize as , e.g. with simple exponential dependencies, Biswas and Hamann demonstrated the computation of to be reducible to order , where is the number of atoms within the interaction range of atom , while irreducible potentials, such as the ATM potential, require a number of operations proportional to .

5.5. Tersoff potentials

Dissatisfied by the poor transferability of the SW potential [235] and inspired by Abell’s analysis of the sensitivity of bond strengths on the local environment [62], Tersoff presented a series of papers [242–246], in which he introduced the concept of bond order to the world of empirical potentials. According to Pauling [247], 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 meant to include nearest neighbors and which then quickly falls off to zero at a distance supposedly less than typical next-nearest neighbor distances:

(68)

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

(69)

where the angular dependence, originally chosen as

(70)

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

Tersoff then constructed a prefactor, , 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

(71)

where , and are element-specific parameters. For large , scales with , which was shown by Brenner [248] to be similar to EAM or TB2M potentials. With these ingredients, a general Tersoff potential reads

(72)

where the (pair) contributions and are strictly repulsive and strictly attractive, respectively. Typically, are exponential functions to reflect the universal equation of state of crystals [62,249].

Using a symmetrized bond-order parameter, , the Tersoff potential can be cast as

(73)

Tersoff chose , as to reproduce the original Morse potential [245]. From the Ducastelle-potential perspective, i.e. when setting 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) 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 in the Tersoff potential, it simply seems wrong that the variable 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 Eq. (53)(53) and (54(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,

(74)

with , where , and are the first, second and third argument of . 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. 2. Both ATM and SW potentials have a many-body contribution of the form , which reduces the many-body contribution to a true three-body contribution, e.g. in Eq. (2)(2) . Carlsson [19] dubbed these types of expressions cluster potentials. Ducastelle and Tersoff introduce a nonlinear mapping for the functional dependency on in , 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(2) ), cluster functionals implicitly include higher-order terms. This can be seen by expressing the many-body contribution as a Taylor series in . This Taylor series is truncated for cluster potentials but contains quadratic , cubic and higher-order terms for cluster functionals. The quadratic term,

(75)

clearly contains a four-body contribution to in Eq. (2)(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(2) ), while avoiding computationally expensive, explicit calculations of higher-order terms.

Table 3. Selected many-body potentials as cast into the generic functional form, Eq. (53)(53) and (54(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. and are characteristic length-scales of repulsion and attraction, respectively. are cutoff functions whose expressions differ between Stillinger-Weber and Tersoff, see Eq. (132). The functions and encode the angular dependence of Stillinger-Weber and Tersoff, respectively. Divergent onsite terms 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 , for modeling hydrocarbon chemistry [10,250]. This potential is now known as the reactive empirical bond-order potential (REBO). Brenner realized that Tersoff’s potential essentially ignores the bond in -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

(76)

where is the coordination of atom obtained with a smooth cutoff function , . The function accounts for the missing -bond and is a bicubic interpolation between values tabulated for specific coordination of the atoms. is another integer variable that is if the bond is conjugated, which is the case if a neighboring atom (of the bond -) has a coordination less than four. increases the neighbor shell, which influences the energy of bond -, making the potential longer ranged than Tersoff’s.

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

(77)

where is again a lookup table and and are coordination numbers of atom , but only counting hydrogens or carbons, respectively, such that . is only nonzero for C-C bonds if a hydrogen is present as a bonding partner, . The lookup tables were fitted to atomization energies of many different hydrocarbon molecules. More information can be found in Refs. [10,250], including nodal values for the lookup table and .

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 and . 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 , depending on the coordination of atom . 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 [251] (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 [252]) 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

(78)

where the is the Fermi-Dirac distribution, are the energy eigenvalues and is the density matrix.² At low electronic temperature, is simply a step function such that the sum in Eq. (79)(79) runs over all states with energy below the Fermi level . The total energy of the system typically includes a pair-wise repulsive contribution in addition to , 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]

(79)

where is orbital on atom . In terms of this local density of states, the band energy becomes

(80)

The moments of the local density of states are given by [253,254],

(81)

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

(82)

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

(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 . 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 [255,256] (that lead to the systematic development of bond-order potentials [257,258]), the maximum entropy method [259–262], or simply approximating the local DOS with a rectangle [202,263]. Let us now assume we have an unary (for simplicity -valent) system (e.g. Cu) with bond integral for the ss-bond. There is no angular dependence and the second moment of the local DOS at atom – see Eq. (83(83) ) – is given by

(84)

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

(85)

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

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 [265] and transition metals [266]. While higher-moment atom-centered techniques have seen some development since [267,268], 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 where the bond energy

(86)

and the promotion energy 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 [269]. Equation (87)(87) for the bond energy is appealing, because it introduces the bond - as the central quantity: It has the appearance of a pair interaction with pair energy ; compare it for example to the attractive part of the Tersoff potential given by Eq. (73)(73) . The term is called the bond order. As we know from Tersoff’s potential, this bond order – or rather the bond order variable , Eq. (72(72) ) – 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 [258].

We can write the density matrix as

(87)

which is easily seen by inserting this expression into Eq. (79)(79) . Note that written this way, the density matrix has a form similar to Eq. (80)(80) and (81(80) ). The onsite terms correspond to the integral over the local density of states up to the Fermi level, which are the number of electrons in that orbital. We now consider the linear combination of orbitals that correspond to bonding () and anti-bonding () states [203]. Then and hence . The bond-order can therefore be written as

(88)

with . This is the difference of electrons in anti-bonding () and bonding () orbitals. We conclude that Pauling’s interpretation [247] of the bond-order was indeed correct.

Exploiting the symmetry of the - bond, the off-diagonal Hamiltonian block can be diagonalized (by rotating into the frame of the local bond -), allowing the pair energy to be written as [270],

(89)

the sum of , and contributions to the bond energy. Pettifor, Aoki, Horsfield and coworkers [257,258,271,272] developed a systematic expansion of the bond-order in the moments computed on atom and . To second order, they find , leading to an expression similar to Tersoff’s total energy (if -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 [273–275], molybdenum [276], gallium-arsenide [277], silicon [278], tungsten [279], iron [280], cadmiun-telleride [281], and aluminum-copper-hydrogen [282–284].

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 [with a single (steep) minimum in 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

(90)

Values for , , are stated in Table 2 for mono-atomic cubic crystals along with other characteristic, dimensionless numbers, which will be introduced further below. Given that is the maximal coordination number, face-centered cubic (fcc) and hexagonal closed packed (hcp), which has the same , , and 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 , where 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 times the energy per atom in the dimer, which is . 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 , despite counteracting effects of other shells. The inequalities

(91a)

(91b)

(91c)

encapsulate this discussion, being the binding energy per atom (in the ground state with nearest-neighbor distance ) and is the volume per atom. Values of for different structures are tabulated in Tab. 2. These equations are expected to hold for non-ionic solids interacting through pair potentials. is the curvature of the pair potential in the minimum, which could be called the dimer spring constant. It takes the value

(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 can be of order so that the sum over shells is not convergent. Summing up Coulomb interactions (in periodically repeated systems) is an entire industry by itself [285,286]. However, as long as crystallographic positions are fixed relative to the unit cell, it is clear that the sum over terms is a geometric factor depending on the lattice times . To keep the treatment analytically amenable, we only consider nearest-neighbor repulsion plus Coulomb interaction and approximate the energy per atom as

(93)

where is the geometric factor for , the Madelung constant. Since 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. (92a)(91a) . Eq. (94)(93) is not necessarily sufficiently accurate to determine what phase an ionic crystal assumes as a function of the parameter [48,150]. To this end, it is necessary to also include dispersive interaction and repulsion between anions plus potentially many-body terms [150], 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 [153],

(94a)

(94b)

(94c)

where follows from Eq. (93)(93) . 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 , the equation of state (EOS) at small temperature follows from , which can be evaluated quite easily through

(95)

when using the volumes per atom 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 Figure 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, , 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 Figure 3 and with parameters listed in Table 1. The pressure and the volume are normalized to the zero-pressure bulk modulus and volume , respectively, which the fits were constrained to reproduce. Black full circes [287] and open circles [288] show reference data for fcc argon, orange plus signs [289] and crosses [290] for fcc copper, red diamonds [295] for carbon in the cubic-diamond structure as well as blue, full [291] and open [292] squares for NaCl in the rock salt structure.

The apparent dimer binding energies deduced from the fits are decreased by a factor of approximately 2.2 for carbon () and by a factor of 4.6 for copper (). 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 .

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

(96)

Here, is the bulk modulus at and , while is the change of the bulk modulus at and constant temperature, which is ideally zero temperature when gauging interaction parameters from lattice sums. In fact, the reference data for carbon [295] in Figure 10 was produced using Eq. (97)(97) 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 parametrizations of potentials make use of such data [296–299].

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(60) ), can be approximated as

(97)

It follows for the equilibrium bond length that

(98)

and for binding energy

(99)

with . Moreover, the generalized stiffness of the potential turns out to be

(100)

where 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 Figure 9—as well as for a few other metals condensing in fcc. These trends are a logarithmic increase of with and an algebraic increase of with a power less than . It is important to note that for we obtain , and Eq. (99)(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 Figure 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

(101)

of nearest-neighbor Tersoffium, Eq. (74)(74) . We directly see that we recover Eq. (98)(98) of Ducastellium for , illustrating the relationship between the bond-order variable and the coordination number, see also Ref. [248]. The bond-order variable typically does not depend on bond-length, but only on crystal structure. We can therefore read-off the equilibrium bond length as,

(102)

Similarly, the energy per bond scales as

(103)

explaining the rough structural trends in Figure 9 for carbon. We can combine Eq. (103)(103) and (104(103) ) to eliminate , yielding

(104)

which is displayed in Figure 9h.

These analytical equations are useful in parameterization of these potentials. Albe [296,297,299] developed a parameterization procedure, based on data of ground-state crystalline structures such as the one shown in Figure 9. The first step is fitting the parameters to the universal energy-bond relation given by Eq. (105)(105) . Albe calls the corresponding plot, such as Figure 9d for copper and 9h for carbon, a Pauling plot. The remaining (pair) parameters , and 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 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, , for which the lattice in the vicinity of the removed atom is kept undistorted, and the true defect energy, , which is obtained after lattice relaxation. The reduction of the defect energy due to relaxation can be deemed irrelevant when studying trends, in contrast to the 10 to 10⁵ 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 – see also Eq. (100(100) ) –, where, however, for pair potentials and for Ducastellium. Removing an atom from an ideal lattice site and placing it hypothetically into another ideal lattice site makes its ditched neighbors lose one bond each, leads to an estimate of, , which becomes

(105)

where the latter approximation applies if . Thus, for a pair potential (), while for , our back-of-the envelope calculation provides a lower, yet good estimate for . Taking into account that the bond lengths of the ditched neighbors is not at the perfect bond length for their new coordination would lead to a relative increase of of typically less than 5% in a Ducastelle potential. Figure 11 confirms the predicted trend that is distinctly less than unity for simple metals, in fact their falls slightly below 0.5, which is the upper bound for approached in the limit , mentioned in the text near Eq. (99(99) ). Group 14 elements have substantially larger than metals, however, and cannot be expected to be similar, whenever bond-angle corrections matter. Only argon out of the elements considered in Figure 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 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 [300]. Vertical defect energies are provided for Al [301], Au [301], Pt [301], Cu [301], Pb [301], Pd [301], Na [302], Ta [301], W [301], Cr [301]. For C, Si and Ge, values were calculated using DFT-LDA. In addition, the vertical defect energy reduction of Ar was deduced from the ATM parameters from Ref. [147].

One reason why defect energies in metals are important is that they are believed to correlate with the melting temperature : Górecki found the vacancy concentration just below in many different metals to be [303] . Estimating the relative number of vacancies with would yield , which roughly overestimates the melting temperatures reported in Figure 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, , the temperature has to be identified at which the Gibb’s free energy of the condensed phase equals that of the gas. At low temperature, can be crudely estimated with the cohesive energy of the crystal and at high temperature with times the entropy of an ideal gas. This gives a quasi-linear relation between and with corrections logarithmic in and . When checking data for metals, we found them to obey a , 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 , when expressing its cohesive energy of [304] in ice Ih per molecule (pm) rather than per atom.

Given the discussions on melting and boiling temperatures, their ratio must be expected to correlate with the exponent . For our representative pair-potential element argon, melting and boiling temperature differ only by 4% at atmospheric pressure. In contrast, , , 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 and should be bad. For reasons of brevity, we do not only abstain from a more detailed analysis of the effects of the 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 .

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

(106)

The first constraint to be discussed is the ratio for fcc crystals interacting with short-range potentials. Using the results compiled in Eq. (35), a first estimate of is obtained using the Cauchy relation and valid for cubic systems interacting through pair potentials. This is already close for to the value of 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 relative to and thus so that the exact value for is increased compared to the nearest-shell approximation. First, the effective spring constant 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 becomes negative. Numerically, we obtain and when using 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 , pushing the nearest neighbors more deeply into the repulsion and leading to a reduction of by and an increase of and thus by the same amount, see also Eq. (39)(39) .

The second constraint to be discussed for a short-range potential is the stability of the bcc phase. In the nearest-neighbor approximation, , which means that the tensor of elastic constants is not positive definite thereby violating the Born mechanical stability criteria [140]. Including the second shell helps stabilizing bcc, but only if the interaction is sufficiently long ranged so that is positive. For a Morse potential, we find to be less than , once and the binding energy per atom for , 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 [158].

Some solid-state physics textbooks argue that a positive next-nearest neighbor coupling stabilizes the simple cubic structure against shear, polonium being the single element in the periodic table adopting it as ground state. However, for to be positive at a next-nearest-neighbor distance of , the Morse exponent has to be as small as , 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(57) ), elastic constants of simple crystals can be found quite easily with the formalism laid out in Sects. 3.5,

(107)

where and are the first and second derivative of evaluated at , respectively, and . In the nearest-neighbor approximation, this expression simplifies to

(108)

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

(109)

in the nearest-neighbor approximation. Thus, if the embedding function has a positive curvature at , as is the case for a Ducastelle potential with , then follows for cubic systems. The symmetrized shear modulus, , where the so-called tetragonal shear modulus is defined as , can be shown to obey in a nearest-shell Ducastelle potential, where for fcc and for bcc [138]. Thus, the most generic EAM or TB2M potential predicts to scale linearly with , just like the vacancy energy. In fact, simple metals obey this trend, as demonstrated in Figure 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 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 for Ar from Ref [305].. Vertical defect energies are provided for Al [301], Au [301], Ni [301], Cu [301], Li [302], K [301], V [301], Nb [301]. The vertical defect energy reduction of Ar was deduced from the ATM parameters from Ref. [147]. For C, Si, Ge and Sr, vertical defect energies were produced using DFT-LDA, with 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” [306]. 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 , it can be seen that the low limit of bring the Poisson’s ratio close to that of a liquid, i.e. to , while the lower limit for is just a little above the value valid for pair potentials . Thus, the smaller , the more “liquid-like” the elastic tensor, which, of course, also implies that the energy to introduce a dislocation, which is linear in , is relatively small. In fact, the elemental solids with a small ratio, but also the ones with a small , are by and large more ductile than the ones with a larger . 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 , , and so on [307,308], 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,308]. 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 (, typical for copper, with 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 graines remained unaltered during the relaxation. Results are shown in Figure 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. [222]. Simulations were run using LAMMPS [380,381]. Colors indicate the orientation of grains, i.e. the angle obtained when averaging over all atoms within a cutoff, where denotes the angle of a bond w.r.t. -axis. Atoms, for which, , are depicted in black.

Figure 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 H has a binding energy of 4.74 eV. The Coulomb interaction in a NaCl molecule with a bond length of Å is about 6.1 eV, after subtracting sodium’s ionization energy, eV, and adding chlorine’s electron affinity, 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 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 [88,91].

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, 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 [309] and Löwdin [310]), (ii) the electronic density (Bader [311] and Hischfeld [312]), (iii) the electrostatic potential (ESP), (iv) displacement-induced changes in polarization (Born [75]), and frequently forgotten (v) experiment [313]. These methods will be sketched next. For more detailed overviews, see Refs. [313–315].

The Mulliken [309] 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 [310] 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 [311] 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 [316]. Hirshfeld [312] 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 [312,317]. In such “iterative stockholder approaches”, it remains challenging to identify good trade-offs between the ESP accuracy and the transferability of the charges [318]. 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 [319]. 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 charges, which can lead to charges exceeding the formal oxidation number. Central ways to deduce partial charges from experiment [313] 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 w.r.t. the electron density disappears, because minimizes . In a coarse-grained description of the electron density in terms of a set of atomic charges, , this minimization principle translates into what Mortier and varying co-authors [320,321] called the electronegativity equalization principle. Assuming that can be expanded into powers of atomic charges, the ground state energy can be written as

(110a)

(110b)

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

The chemical meaning of the parameters , , and is best ascertained in limiting cases. Consider atom to be neutral and isolated. If one electron is added to it so that , the energy decreases by its electron affinity , while it increases by its ionization energy when one electron is taken away from it. Thus,

(111a)

(111b)

so that the so-called chemical hardness [322] turns out to be and the electronegativity [323,324] . Meanwhile, can be associated with the Coulomb potential given the charges and of two distant atoms and :

(112)

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

Rappé and Goddard [325] 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

(113)

where the charge densities are those associated with -type Slater orbitals. This way, the Hessian, i.e. the matrix formed by 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 elementary charges when truncating the expansion of Eq. (111) after the quadratic term [326].

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 [326], DFT having related issues, QEq makes the polarizability of polymers grow superlinearly rather than linearly in the chain length [327,328], and, similarly, solids adopt the dielectric response function of a metal [329], 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 [330].

One fundamental problem of Eq. (111) 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 [331]. This property should be reflected in any coarse-grained approximation to DFT, which is built on the Hohenberg-Kohn theorem. However, Eq. (111) 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 [329]. 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 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 [332]. For nanostructured electrodes, the finite density of states at the Fermi level leads to an additional contribution to the capacity of the device [333], which can be modeled by replacing the quadratic term in Eq. (111) with an appropriate nonlinear function [334].

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 [335], 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. [327] introduced atom-atom charge transfer (AACT) variables so that the charge of an atom is given by

(114)

where is the charge transferred from atom to atom . In the AACT method, the bond hardness, , replace the atomic hardness used in regular QEq so that a charge transfer costs a potential energy of . 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, , turn out to be barely exceed unity [336]. Moreover, AACT suppresses the observed chain-length dependent polarizability of short oligomers [328,329] and, similarly, produces a zero screening length in dielectrics [329]. 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. [336] 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”, , formerly known as AACT variable, is split between two atoms. Thus, their model reads

(115)

where the entering , i.e. Eq. (111), are defined as in Eq. (115)(115) . Including finite bond hardness suppresses the superlinear scaling of molecules with linear size [337], as is shown in Figure 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 [326]. 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 [154], 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) [338]. 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 matrix must be defined, where is the number of charges, while the number of split charges is only of order . 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 diverges at large . This limitation can be overcome by introducing oxidation numbers so that an atomic charge is given by

(116)

with in a so-called redox-SQE formalism [157,339]. 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 . The relation between point-dipole polarizability and is text-book material and reflected in the Clausius-Mossotti (CM) relation, Eq. (42(42) ), while that of a bond-polarizable model without point dipole polarizability was worked out by Nistor and Müser [329]. 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 [340]. It starts from the insight that the full electrostatic field of a dipole contains a contribution 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 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 (which is assumed to be constant) minus , where is the density of point dipoles. For a rocksalt lattice, or, when discretizing a homoegeneous material into cubes of linear size and arranging the local charge transfer variables, or, split charges in a vector , the dipoles obey associated with the charge-transfer, and the point dipoles obey

(117a)

(117b)

respectively. The x-component of q is the split charge donated would be the split charge donated from an atom or grid point to its right neighbor, and so on, while is the corresponding bond hardness. The total polarizability is nothing but the total dipole density , while is defined through within linear response. After some minor algebra

(118)

can be obtained, which contains the CM relation in the limit of as well as valid for bond-polarizable models ignoring point dipoles [329]. Once bond polarizability is ignored in a charge-transfer model, 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 . A wave-number dependent dielectric constant and thus a dielectric smearing or decay length ensues. For a pure bond-polarizable model and the geometry considered above, [336]. One of the consequences of this result is that ideal mirror charges do not accurately represent the dielectric response of real metals [332].

The potential importance of modeling dielectric response correctly is also demonstrated in Figure 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. Note that Figure 15 is a proof-of-principle demonstrator for the redox-SQE method [339]. 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 [341].

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).

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 methods.³ 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 [342]. 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 [343–346].

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 [238]: “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 [242] turned out to have bcc and not diamond cubic as the ground state [347]. Yet, potentials constructed on carefully curated databases give excellent transferability even to obscure crystal phases, see e.g. Ref. [348] for a current example, again for elemental silicon. This highlights the specific importance of database quality for MLPs.

Most MLPs consists 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. [29,31,32,349,350]. 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. [350,351], 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,11,250]. Given a bond -, Brenner used the coordination numbers , and a bond conjugation parameter 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, where 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 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 -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)(2) . An objective descriptor for two-body interaction is simply the bond distance . For a three-body term, we can use the triplet of distances as a descriptor. While this descriptor satisfies objectivity, it still has permutational symmetries with respect to the exchange of the positions of like atoms and . Bartók and Csányi [349] suggested the use of the descriptor

(119)

which is equivalent to but symmetric with respect to the permutation .

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 triplet, can also be written in terms of the descriptor . 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. (120)(120) .

The specific descriptor of Eq. (120)(120) was used in the Gaussian approximation potential (GAP) [352] for amorphous carbon [353], later extended for improved accuracy of crystalline phases [354]. 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 as a descriptor. The second step involved correcting this potential using a three-body term employing the descriptor given by Eq. (120)(120) . 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(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 [355,356]. The atom-centered fingerprints consist of radial symmetry functions

(120)

These functions are Gaussians centered at a pair distance , which are forced smoothly to zero at a cutoff distance via . They are similar to coordination numbers, but put specific weight on atoms at the distance . In addition to these radial functions, there are angular functions

(121)

We have added the superscripts and to emphasize, that those functions are not evaluated for a single parameter vector and cutoff function but for a set of these parameters. The fingerprint vector is then constructed by combinations of and 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 and drop the double notation , , implying that the index also refers to the type of symmetry function and not just its parameters.

To emphasize connection to classical potentials, the look like the onsite terms occurring in our general functional form, Eq. (54(54) ), while the look like offsite terms. Indeed, if we replace the Gaussian in Eq. (121)(121) with an exponential, the expression becomes identical to the local density of Ducastellium, see Eq. (59b)(59b) . Similarly, 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 [357]. It can be regarded as a generalization of the bond-orientational order parameter of Steinhardt, Nelson and Ronchetti [358]. The idea is to map the atomic configurations onto an atomic density field , which is given by atom-centered functions,

(122)

where is broadened to a Gaussian with a characteristic length . Note that 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 for each atom ,

(123)

with local angles and of the vector . The functions are orthogonal radial basis functions. The power-spectrum of the angular coefficients ,

(124)

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

SOAP taken to higher body order – and with as in Eq. (123(123) ) – forms the basis of the atomic cluster expansion (ACE) [359,360]. We recommend study of Ref. [359] 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 [361].

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 , such that the total energy is simply . 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 [362], which uses the fingerprints introduced in the original GAP [352]. A more flexible regression strategy is the use of Gaussian process regression [252], which lends its name to the GAP potential [352]. 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 [363]. The final expression for the per atom energies takes the simple form,

(125)

where the sum over runs over all elements of the training database. The fingerprints of the current environment of atom are compared to all fingerprints in the training database. The function 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 [252]). From Eq. (126)(126) 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 is a squared exponential (Gaussian) [252,353], which containsa set of distance scales as a hyper parameter,

(126)

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

An interesting use of Gaussian process regression is the on-the-fly parameterization [365] of a force field based on electronic structure calculations carried out on snapshots throughout a molecular dynamics trajectory. Li, Kermode and de Vita [366] 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 [365], e.g. of fracture [367,368], 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 [359,369] 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 [369]. Fully nonlinear regression can be achieved with neural networks. Behler and Parrinello [355] employed the specific form

(127)

for regression in their neural network potential. Here, and are sigmoidal “activation functions” and and are weights. The weights are numbers that are adjusted in the process of fitting the neural network. The sum over runs over all entries of the feature vector of atom i. This network has a single “hidden layer” with activation function and the index runs over the outputs of this layer. The original fit for silicon by Behler and Parrinello employed three nodes in the hidden layer and entries in the fingerprint vector, amounting to a total of weights. Equation (128)(128) resembles the many-body contribution of classical potentials for open-shell systems, see Eq. (53)(53) .

9. Limiting the interaction range

Chemistry is local for the most part [370]. 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 [371–373]. 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 [225,374] 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 . He continued that “this surprising result gives some support to the procedure we have used here in ignoring all but first-neighbor interactions” [225] (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 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 [375,376] and thereby make it on the cover of contemporary text books [377], or, nano-tube-based water pumps work endlessly without external driving due to improper cut-offs [378], the report of which [379] 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 and for the efficient parallelization on parallel, high-performance computing systems using domain-decomposition [380–386]. Even long-range Coulomb potentials can be computed with a numerical effort scaling linearly or quasi-linearly with [387–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 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 -function singularity in the interaction force at , 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 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 passes through , 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 such that all its derivatives up to ‘th order approach zero continuously at and remain zero for :

(128)

Here, is the ‘th order Taylor expansion of about .

The effect of cut-and-shift procedures on and on the equation of state is shown exemplarily for a Morse potential (, , ) in Figure 16a. It becomes evident that “large” lead to a poor representation of the binding energy unless is large, which, however costs computing time. However, serious artifacts can occur when potentials are truncated but also when they are merely shifted with . The EOS of an athermal crystal is discontinuous even for as revealed in Figure 16b. Thermal fluctuations smear out and potentially eliminate this discontinuity. Nonetheless, caution remains advisable. In principle, a ‘th-order shifted potential can have a discontinuity in 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 and a cut-off distance of . The cut-and-shift procedures include (blue), (red), and (green). Right panel: energy per atom of an ideal fcc crystal as a function of the nearest-neighbor distance using the shifted potential (blue), shifted force (red), and shifted curvature (green) for the Morse potential with and a cutoff radius . Dots and squares denote reduced volumes at which a neighbor-shell radius is equal to .

An alternative to cut-and-shift potentials is the use of cut-off or smoothing functions , which are multiplied with the true pair potential to yield

(129)

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

(130)

with . 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 with which derivatives of the potential disappears at . However, an additional risk of smoothing functions is that using “large” does not guarantee artifacts to be suppressed in particular when is close to . Prominent smoothing functions are [216,217,224,227,242,393]:

(131a)

(131b)

(131c)

(131d)

(131e)

(131f)

where integer Mishin exponents are commonly selected as . In Eq. (132), and specify which highest ‘th-order derivative of the smoothed function approaches zero at and , respectively. Moreover, means that while , while omitting means that . Recently, it was found that cutting with polynomials constructed such that 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 than at .

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 contribution (which could be seen as a kind of mean-field correction), and to ignore any contribution at non-zero wave vector . Some authors find support for the correctness of Wolf’s optimistic assessment [285,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 in silica melts [394]. However, problems arise when structural heterogeneities extend on length scales [394]. They induce undulations in chemical potentials with wave vectors of order , which favor a charge separation. Coulomb interactions suppress that separation with terms scaling as . 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 [157], 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 [285,394], can be made for homogeneous systems, but steer clear of it otherwise, in particular in any biophysical context [286].

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 (Figure 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 Figure 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 Figure 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, respectively [222]. The MEAM cutoff function given by Eq. (132b)(131b) was applied to the repulsive potential and the density. and are the cohesive energy and the nearest-neighbor distance for an fcc crystal described by Ducastelle potential, respectively.

9.2. Screening

In order to reflect screening through screening rather than through cutting, Baskes et al. [225,374] proposed to replace cutoff functions with environment-dependent screening functions, , which themselves are products of three-body functions ,

(132)

The are constructed such that when atom screens the interaction between and completely, while when atom does not screen the interaction at all. Thus, the screening function acts like a topological cutoff function, . Baskes proposed different functional forms for . The refined expression is [374]

(133)

which drops from unity to zero between and . The quantity characterizes the geometry of the ellipse that passes through the three atoms --. It is given by

(134)

the square of the ratio of the lengths of the two half-axes of that ellipse. Here is the distance between atoms - normalized by the length of the central - bond.

Figure 18 illustrates this concept. The bond is unscreened () if there is no third atom in its vicinity (Figure 18a). As a third atom moves into the shaded ellipsis (Figure 18b – the geometry of the ellipsis is described by the values of and ) the bond becomes screened (, Figure 18c). It becomes clear, that if we simply rescale all coordinates (as in the homogeneous, volumetric dissociation of a crystal), remains constant for each bond and an energy would follow the true nearest neighbor curve, see the black line in Figure 17(a). The screening function 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. (a) Initially, the thick, black bond is unscreened, but as (b) the atom enters the region of influence (light gray ellipses) the bond weakens. (c) The atom has moved into the close vicinity of the bond (dark gray region), effectively disabling it while creating two new bonds (red lines). The Baskes screening functions [374] are defined the aspect ratio of the light gray region () and the dark gray (), defining a measure of bond screen that is independent of the absolute lengths of the bonds in the system.

The Baskes screening functions [374] were applied to empirical bond-order potentials independently by Pastewka et al. [25,376,401] and Kumagai et al. [402]. Both groups emphasized that the screened potentials significantly improved the properties of amorphous carbon modeled with REBO2 or Tersoff-type potentials. In addition, the screening functions served to overcome the issue with dissociation of a bond under external stress discussed in Refs. [25,375,376,401]. This enabled modeling of fracture in crystalline and amorphous carbon systems [403]. Perriot et al. [404] presented a slightly different screening concept requiring the REBO potential to be refitted.

Screening functions can be rationalized as originating from nonorthogonality in a tight-binding framework. This nonorthogonality leads to an environment-dependence of the bond-integrals, when the nonorthogonal tight-binding is “coarse-grained” to an equivalent orthogonal tight-binding model. Nguyen-Manh, Pettifor and Vitek [405,406] showed, that the theory of the bond-order expansion, briefly touched upon in Sect. 5.8, can be used to derive screening functions from a nonorthogonal tight-binding model. This first-principles construction lends additional support to Baskes’ screening functions and other screening approaches, such as empirical environment-dependence introduced in the context of orthogonal tight-binding shortly after Baskes work [407–410].

10. Summary and perspectives

Eugene Wigner allegedly said: “It is nice to know that the computer understands the problem. But I would like to understand it too.” One main motivation for writing this review was to assist people with similar ambitions as Wigner. To this end, we summarized our understanding of what properties in condensed-matter systems can be induced by the functional form of the potentials used for their description. In our endeavor, we felt compelled to create much own data and new figures with the purpose to create insight and to convey trends and differences between potential classes rather than to produce numbers for a specific system. When doing so, we did our best to embed anything written into a historical context, which is summed up in Figure 19.

Figure 19. Selected highlights of the development of interaction potentials.

Yet many times, we could not find appropriate references or quotes, which we are certain do exist. As one of numerous examples, we found many papers computing shear and bulk moduli of metals or vacancy-defect and cohesive energies, but always missed the argument why their respective ratios are correlated and how they relate to the ratio of melting and boiling temperature. We expect our discussion to have satisfied Wigner’s desire for understanding the correlation between these ratios. Despite certainly having missed well known studies, we did find some old works, which may have been underappreciated, such as Slater’s paper on the interaction between helium atoms [82]. As mentioned earlier, Slater derived the exponential repulsion between atoms with closed valence shell, which promoted Born and Mayer [83] as well as Buckingham [7] to use this or slightly modified forms for the repulsion in the potentials now carrying their names. However, Slater also derived the dispersive coefficient for helium to within 15% accuracy, two years before London [56] generalized the results to other closed-valence shells.

Despite the length of this article, we could only scratch the surface of the large field of interatomic potentials. Many central aspects were not touched upon. Most importantly, we barely discussed how to adjust parameters, in particular the pros and cons to fit to experimental or to in-silico data. There are good reasons to follow the main-stream opinion that the quality of a potential increases the less empirical the data on which the potential is parameterized. Computer-generated reference data is much more versatile than that provided from experiments. Forces on individual atoms can be used for characteristic bonding situations or rare but important configurations like a transition state occurring during a chemical reaction or a collective phase transformation [411]. Moreover, in-silico data does not contain quantum effects, which frequently need to be accounted for when comparing computer-generated data to experiments. Describing how to do that properly would have required us to outline how to approximately subtract the nuclear quantum effects from experimental data or how to incorporate quantum fluctuations into the simulations, e.g. through path-integral techniques, or, by encoding their effect into effective temperature-dependent, many-body potentials, which would have been beyond the scope of this review.

Thus, there is scarcely any argument to gauge parameters on experiments, if there was not the small but important detail that experimental data is by and large more accurate than density-functional theory, which cannot be deemed exact, as long as the exact functional has not been identified. It could be argued that we base one theory or potentials with uncontrolled approximations on another one with uncontrolled approximations, which has trouble to predict two dislike molecules or clusters separated by a large distance to each acquire an integer charge [412,413]. When dismissing empiricism as fundamentally problematic, one may also keep in mind that one of the greatest theoretical achievements in chemistry, arguably in all of science, was the construction of the periodic table by Mendeleev. He even predicted the existence of unknown elements including some of their physical and chemical properties with an accuracy that people using potentials or even DFT might have a hard time to match if they did not know what they had to predict, or, rather postdict. Moreover, the amount of data that Mendeleev could build on was noticeably less than what is required in machine learning.

The potentials discussed in this review pertain mostly to situations, in which bonds can be clearly classified as dispersive, metallic, covalent, or ionic. For situations, where this simple categorization cannot be made, different potential classes are combined in a mix-and-match fashion into compound potentials. Prominent examples are the adaptive intermolecular REBO (AIREBO) [251] (combining Brenner’s potential with nonbonded Lennard-Jones interaction), the Streitz-Mintmire potential [414] (combining EAM with charge transfer), the charge-optimized many-body potential (COMB) [415–417] (combining Tersoff’s potential with charge transfer), a merger of REBO with split-charge equilibration [330,418], as well as early combinations of Keating-type with charge-transfer potentials potentials [419,420]. Of course, the widely-used ReaxFF potential [12,421], which merges a bond-order approach (different in nature than the approaches discussed in this review) with non-bonded interactions and charge transfer, must also be mentioned.

While compound approaches can be extremely powerful, many of them simply add different energy terms. This can be problematic even for seemingly simple alloys or intermetallics formed by elements of large electronegativity difference. Put simply, negatively charged atoms grow in size while positively charged atoms shrink. This symmetry breaking between negative and positive charge is not reflected when simply adding charge equilibration to a (post) Ducastelle potential. Yet, it is supposedly responsible for why the negatively charged atoms in intermetallics have the tendency to close pack while positive atoms occupy interstitial positions, as it happens, for example, for AlAu, also called the purple plague: Au atoms form an fcc lattice while Al atoms assume interstitial positions. Although promising steps toward true compound potentials have been taken [422], e.g. by augmenting or reducing the valence density of a neutral atom with a term proportional to its partial charge, systematically merged potential remain a dream.

Novel paths that are taken with machine learning potentials seem extremely promising. However, a puzzling question is why machine-learned potentials outperform parameterized potentials. The claim that they are parameter free or free of functional constraints is not entirely justified. Many of the local descriptors are suspiciously close to what is used in potentials, as indeed they are often “physics-inspired” [350]. However, the big advantage of MLPs is that they do not make strong assumptions like pair-wise additive repulsion, which might be one of the most important sources of error in classical interatomic potentials.

A show-stopping problem central to all potentials is the curse of dimensionality. Fitting multi-species (or alloy) potential requires a number of pair-parameters that scales asymptotically as with the number of atomic species . The scaling becomes even less favorable if we need specific parameters for triplets (), quadruplets () and so on, quickly becoming intractable for a large number of species. The compression of chemical fingerprints has recently been proposed to circumvent the curse of dimensionality for MLPs [423]. Using explicit functional forms, it can be possible to circumnavigate the curse of dimensionality with combining rules. However, they are only available for few interactions types and may be plagued with poor transferability.

As a final note, we would like to point out that despite the fact that (with the exception of bare Coulomb interaction) all potentials discussed here are local, chemistry can be quite non-local. By non-locality we do not mean the range of the bare interaction, such as the range of the bond-integrals in a tight-binding formulation. We mean the non-locality intrinsic in the diagonalization of the quantum mechanical Hamiltonian. In hydrocarbon chemistry, the non-locality manifests itself for example in bond conjugation and in metals through an algebraic decay of the density matrix [424], while in group 15–17 in the periodic table, it is reflected in the Peierls deformation causing elemental crystals to reduce from the simple cubic to less symmetric structures [425]. As another example, carbon chains – also called carbynes – can exist in a polyynic form of alternating single and triple bonds or a cumulenic form of repeating double bonds [426–428]. Which form is chosen depends on whether the chain is odd or even numbered and how it is terminated. This crucially affects how they interact with their environment, for example with oxygen [429,430]. Such non-local effects even manifest in bulk materials: Force-locality tests on amorphous carbon by Deringer and Csányi showed that chemistry in low-density, graphite-like amorphous carbon is much longer ranged than in denser more diamond-like carbon [353]. Approaches for incorporating true quantum non-locality into potentials currently do not appear to exist. Modeling it appears to require new classes of potential, e.g. the ability of an EAM or MEAM potential to make atoms adjust their donating charge density in response to the environment in a fashion that allows for multistability.

We hope that this review was successful in highlighting the incredible achievements throughout the last century in understanding the bonding of matter, and molding these insights into simple analytical expressions. The wide availability of high-accuracy electronic structure calculations and advances in statistical modeling have moved the field into exciting new directions. We would also like to add that the wide availability of present-day interatomic potentials in the form of open-source software, ideally embedded in a standard database [13] or a standard code [380,381], is accelerating quick adoption of potentials into practice — not to mention the savings in students’ lifetimes, by not having to dissect which of the parameters just manually copied from printed publication XYZ is missing a in print. (Yes, we are thinking about our own PhD theses.) Of course, significant challenges remain, both for traditional fixed-form as well as machine-learned interatomic potentials, of which we believe the curse of dimensionality and the coupling of electron transfer and Coulomb interaction to the electronic bond as most crucial [431].

Nomenclature/Notation

Table

Abbreviations
ACE	atomic cluster expansion
ATM	Axilrod-Teller-Muto interaction potential
DFT	density-functional theory
EAM	embedded-atom method
EOS	equation of state
GAP	Gaussian approximation potential
LJ	Lennard-Jones
MEAM	modified EAM
ML	machine learning
MLP	machine-learned potential
QEq	charge equilibration
REBO	reactive empirical bond-order potential
REBO2	second generation REBO
SQE	split-charge equilibration
SW	Stillinger-Weber
TB, TBnM	tight-binding, TB n-th order moment expansion

Symbols

, polarizability in SI and atomic units

Madelung constant

Kronecker delta

Lennard-Jones energy parameter

vacuum permittivity

dielectric constant

element of the Eulerian strain tensor

element of the Lagrangian strain tensor

charge density, number density

length scale parameter

element of the Cauchy stress tensor

electron affinity

-th order virial coefficient

element of elastic tensor

element of elastic tensor in Voigt notation

dispersion coefficient of order

Hamilton operator

Hamiltonian integral between orbital on atom and orbital on atom

ionization energy

bulk modulus

particle number

charge of atom

square of the distance between atoms and

interaction energy

dimer/molecular binding energy

potential energy per atom

equilibrium potential energy per atom (cohesive energy)

potential energy per bond

equilibrium potential energy per bond

single-body interaction energy

pair and triplet interaction energy

volume

coordination number

number of atoms in ‘th nearest-neighbor shell

distance between an atom with a -nearest neighbor

equilibrium bond length

Bohr radius

cutoff function

thermal energy

mass

density of states

local density of states of orbital

momentum operator

dipole moment of species and its magnitude

pressure

descriptor of the environment of atom

bond charge donated from atom to atom

, (pair) distance

equilibrium distance in a diatomic molecule

cutoff radius

, second- and fourth-rank shell tensor for the ‘th nearest-neighbor shell

volume per atom

Acknowledgments

We thank Gábor Csányi, Volker L. Deringer, Christian Elsässer, Peter Gumbsch, Judith Harrison, James Kermode, Pekka Koskinen, Gianpietro Moras, Michael Moseler, Matous Mrovec, Toon Verstraelen and Michael Walter for many enlightening discussions over the years. We are further indebted to James Kermode for comments on the machine learning section of the manuscript as well as Joshua Weißenfels and Jan Grießer for proofreading and commenting on the full manuscript. We used gpaw [431] for all DFT calculations shown here that are not obviously taken from third sources.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1. A large set of these tabulated potentials files can be found in the NIST interatomic potentials repository at https://www.ctcms.nist.gov/potentials/.

2. Functions of the Hamitonian are defined through their action on eigenfunctions, , or through a power series expansion of .

3. See https://libatoms.github.io/GAP/data.html for links to high-quality ab-initio databases for C, Si, and other materials.