Interatomic potentials are mathematical functions to calculate the potential energy of a system of atoms with given positions in space.M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Oxford University Press, Oxford, England, 1989.Daan Frenkel and Berend Smit. Understanding molecular simulation: from algorithms to applications. Academic Press, San Diego, second edition, 2002.R. Lesar. Introduction to Computational Materials Science. Cambridge University Press, 2013. Interatomic potentials are widely used as the physical basis of molecular mechanics and molecular dynamics simulations in computational chemistry, computational physics and computational materials science to explain and predict materials properties. Examples of quantitative properties and qualitative phenomena that are explored with interatomic potentials include lattice parameters, surface energies, interfacial energies, adsorption, cohesion, thermal expansion, and elastic and plastic material behavior, as well as chemical reactions.N. W. Ashcroft and N. D. Mermin. Solid State Physics.Saunders College, Philadelphia, 1976.Charles Kittel. Introduction to Solid State Physics. John Wiley & Sons, New York, third edition, 1968.
Functional form Interatomic potentials can be written as a series expansion of functional terms that depend on the position of one, two, three, etc. atoms at a time. Then the total potential of the system can be written as
Here is the one-body term, the two-body term, the three body term, the number of atoms in the system, the position of atom , etc. , and are indices that loop over atom positions.
Note that in case the pair potential is given per atom pair, in the two-body term the potential should be multiplied by 1/2 as otherwise each bond is counted twice, and similarly the three-body term by 1/6. Alternatively, the summation of the pair term can be restricted to cases and similarly for the three-body term , if the potential form is such that it is symmetric with respect to exchange of the and indices (this may not be the case for potentials for multielemental systems).
The one-body term is only meaningful if the atoms are in an external field (e.g. an electric field). In the absence of external fields, the potential should not depend on the absolute position of atoms, but only on the relative positions. This means that the functional form can be rewritten as a function of interatomic distances and angles between the bonds (vectors to neighbours) . Then, in the absence of external forces, the general form becomes
In the three-body term the interatomic distance is not needed since the three terms are sufficient to give the relative positions of three atoms in three-dimensional space. Any terms of order higher than 2 are also called many-body potentials. In some interatomic potentials the many-body interactions are embedded into the terms of a pair potential (see discussion on EAM-like and bond order potentials below).
In principle the sums in the expressions run over all atoms. However, if the range of the interatomic potential is finite, i.e. the potentials above some cutoff distance , the summing can be restricted to atoms within the cutoff distance of each other. By also using a cellular method for finding the neighbours, the MD algorithm can be an O(N) algorithm. Potentials with an infinite range can be summed up efficiently by Ewald summation and its further developments.
Force calculation
The forces acting between atoms can be obtained by differentiation of the total energy with respect to atom positions. That is, to get the force on atom one should take the three-dimensional derivative (gradient) of the potential with respect to the position of atom :
For two-body potentials this gradient reduces, thanks to the symmetry with respect to in the potential form, to straightforward differentiation with respect to the interatomic distances . However, for many-body potentials (three-body, four-body, etc.) the differentiation becomes considerably more complex
since the potential may not be any longer symmetric with respect to exchange. In other words, also the energy of atoms that are not direct neighbours of can depend on the position because of angular and other many-body terms, and hence contribute to the gradient .
Classes of interatomic potentials Interatomic potentials come in many different varieties, with different physical motivations. Even for single well-known elements such as silicon, a wide variety of potentials quite different in functional form and motivation have been developed. The true interatomic interactions are quantum mechanical in nature, and there is no known way in which the true interactions described by the Schrödinger equation or Dirac equation for all electrons and nuclei could be cast into an analytical functional form. Hence all analytical interatomic potentials are by necessity approximations.
Over time interatomic potentials have largely grown more complex and more accurate, although this is not strictly true. This has included both increased descriptions of physics, as well as added parameters. Until recently, all interatomic potentials could be described as "parametric", having been developed and optimized with a fixed number of (physical) terms and parameters. New research focuses instead on non-parametric potentials which can be systematically improvable by using complex local atomic neighbor descriptors and separate mappings to predict system properties, such that the total number of terms and parameters are flexible. These non-parametric models can be significantly more accurate, but since they are not tied to physical forms and parameters, there are many potential issues surrounding extrapolation and uncertainties.
Parametric potentials Pair potentials The arguably simplest widely used interatomic interaction model is the Lennard-Jones potential .
where is the depth of the potential well and is the distance at which the potential crosses zero. The attractive term proportional to in the potential comes from the scaling of van der Waals forces, while the repulsive term is much more approximate (conveniently the square of the attractive term). On its own, this potential is quantitatively accurate only for noble gases and has been extensively studied in the past decades, but is also widely used for qualitative studies and in systems where dipole interactions are significant, particularly in chemistry force fields to describe intermolecular interactions - especially in fluids.
Another simple and widely used pair potential is the Morse potential, which consists simply of a sum of two exponentials.
Here is the equilibrium bond energy and the bond distance. The Morse potential has been applied to studies of molecular vibrations and solids, and also inspired the functional form of more accurate potentials such as the bond-order potentials.
Ionic materials are often described by a sum of a short-range repulsive term, such as the Buckingham pair potential, and a long-range Coulomb potential giving the ionic interactions between the ions forming the material. The short-range term for ionic materials can also be of many-body character.
Pair potentials have some inherent limitations, such as the inability to describe all 3 elastic constants of cubic metals or correctly describe both cohesive energy and vacancy formation energy. Therefore, quantitative molecular dynamics simulations are carried out with various of many-body potentials.
Repulsive potentials
For very short interatomic separations, important in radiation material science, the interactions can be described quite accurately with screened Coulomb potentials which have the general form
Here, when . and are the charges of the interacting nuclei, and is the so-called screening parameter. A widely used popular screening function is the "Universal ZBL" one.J. F. Ziegler, J. P. Biersack, and U. Littmark. The Stopping and Range of Ions in Matter. Pergamon, New York, 1985. and more accurate ones can be obtained from all-electron quantum chemistry calculations In a comparative study of several quantum chemistry methods, it was shown that pair-specific "NLH" repulsive potentials with a simple three-exponential screening function are accurate to within ~2% above 30 eV, while the universal ZBL potential differs by ~5%–10% from the quantum chemical calculations above 100 eV. In binary collision approximation simulations this kind of potential can be used to describe the nuclear stopping power.
Many-body potentials The Stillinger-Weber potential is a potential that has a two-body and three-body terms of the standard form
where the three-body term describes how the potential energy changes with bond bending. It was originally developed for pure Si, but has been extended to many other elements and compounds and also formed the basis for other Si potentials.
Metals are very commonly described with what can be called "EAM-like" potentials, i.e. potentials that share the same functional form as the embedded atom model. In these potentials, the total potential energy is written
where is a so-called embedding function (not to be confused with the force ) that is a function of the sum of the so-called electron density . is a pair potential that usually is purely repulsive. In the original formulation the electron density function was obtained from true atomic electron densities, and the embedding function was motivated from density-functional theory as the energy needed to 'embed' an atom into the electron density. However, many other potentials used for metals share the same functional form but motivate the terms differently, e.g. based on tight-binding theory or other motivations.
EAM-like potentials are usually implemented as numerical tables. A collection of tables is available at the interatomic potential repository at NIST (see here).
Covalently bonded materials are often described by bond order potentials, sometimes also called Tersoff-like or Brenner-like potentials.
These have in general a form that resembles a pair potential:
where the repulsive and attractive part are simple exponential functions similar to those in the Morse potential. However, the strength is modified by the environment of the atom via the term. If implemented without an explicit angular dependence, these potentials can be shown to be mathematically equivalent to some varieties of EAM-like potentials. Thanks to this equivalence, the bond-order potential formalism has been implemented also for many metal-covalent mixed materials.
EAM potentialshave also been extended to describe covalent bonding by adding angular-dependent terms to the electron density function , in what is called the modified embedded atom method (MEAM).
Force fields
A force field is the collection of parameters to describe the physical interactions between atoms or physical units (up to ~108) using a given energy expression. The term force field characterizes the collection of parameters for a given interatomic potential (energy function) and is often used within the computational chemistry community. The force field parameters make the difference between good and poor models. Force fields are used for the simulation of metals, ceramics, molecules, chemistry, and biological systems, covering the entire periodic table and multiphase materials. Today's performance is among the best for solid-state materials, molecular fluids, and for biomacromolecules, whereby biomacromolecules were the primary focus of force fields from the 1970s to the early 2000s. Force fields range from relatively simple and interpretable fixed-bond models (e.g. Interface force field, CHARMM, and COMPASS) to explicitly reactive models with many adjustable fit parameters (e.g. ReaxFF) and machine learning models.
Non-parametric potentials