Using topology for understanding your computational results

Introduction

In the past two decades, the scientific community has celebrated major anniversaries in physics and chemistry, including Dalton’s atomic theory, ¹ and Gilbert Newton Lewis’s electron pair concept. ² These milestones highlight how chemistry has evolved from an independent branch of knowledge to a fully connected discipline seemingly rooted in quantum theory (QT). Similarly, the advent of computational methods, derived from the theoretical advances of von Neumann and Turing and the development of high performance computing hardware, has made it possible to find almost exact solutions to the Schrödinger equation for small systems. However, Van Vleck’s prediction of the wave function information problem has come true, demonstrating an exponential growth of computational complexity as the system size grows, which is beginning to be mitigated thanks to the explosion of machine learning solutions. ³

Beyond computational complexity, it has become clear that limitations arise also regarding interpretations. QT does not seamlessly translate into the pre-quantum language of chemistry. While QT provides a well-defined algebraic framework, chemists continue to rely on empirical models that remain closer to Dalton’s or Lewis’ paradigms. Atoms, bonds, and functional groups are central to chemical reasoning, yet these concepts lack a direct representation in the Hilbert space where quantum state vectors reside. This has led Paul Popelier and other theorists to point out that chemistry operates at an emergent level that is not fully reducible to quantum mechanical descriptions, raising epistemological challenges examined by philosophers of science ^{4 , 5 , 6} that are similar to those encountered in other complex systems.

To date, chemistry has extensively incorporated QT, with quantum mechanical models guiding chemical intuition. Paradoxically, however, increasing computational accuracy leads to a loss of unambiguous chemical concepts, a phenomenon noted by R.S. Mulliken as early as 1965. ⁷ Some theoretical chemists have attempted to reconcile these two paradigms by developing frameworks in which the standard chemical concepts emerge naturally within a quantum mechanical context.

One significant breakthrough in this endeavor was Density Functional Theory (DFT), established by Hohenberg and Kohn. ⁸ Their foundational theorem states that the electron density ρ determines the full quantum state of a system in its ground state. Later, Riess and Münch in 1981 ⁹ and Mezey’s holographic theorem in 1999 ¹⁰ strengthened this idea, demonstrating that local knowledge of the electron density suffices to reconstruct the entire wavefunction of a system. This insight supports the use of reduced density matrices to condense quantum information into a chemically interpretable form, but most contemporary chemistry curricula and textbooks use the orbital paradigm instead. Historically, the preference for orbital-based models in chemistry has been clearly driven by their success in explaining experimental data. However, alternative approaches, such as those championed by Bader, Henneker, and Cade, ¹¹ emphasize the direct analysis of electron density distributions rather than arbitrary orbital constructs. These perspectives challenge the conventional orbital view by focusing on physically observable quantities rather than model-dependent representations, and have led to the development of the topological approach in chemistry see, e.g. Ref. 12] for a more comprehensive approach, with interesting consequences on the philosophy of chemistry. ¹³

Extracting chemical insights from the electron density has been a long-standing challenge. Early approaches, such as Berlin’s force analysis, ¹⁴ classified electron redistribution into bonding and antibonding regions. These methods illustrated charge accumulation patterns but lacked generalizability. Density difference maps, which compare atomic (or ionic) and molecular/crystalline electron densities, emerged as a widely used technique. However, they face difficulties, such as selecting an appropriate reference density. For example, in homonuclear diatomic molecules, differing choices of atomic references yield contradictory bonding interpretations. Moreover, promolecular densities violate Pauli’s exclusion principle, rendering them inconsistent with quantum mechanics.

After decades of development, the so-called topological approach (also known as quantum chemical topology, see Ref. 15] for a recent review in the context of electron densities) has reached a considerable level of maturity, leading to interesting philosophical discussions about the ontological status of chemistry itself and its relations to physics. Be that as it may, the approach can now be safely used by non-specialists. Although it has branched out into many different flavors over the years, the broader chemical community has come to a sort of consensus, which may or may not be shared by modern theoretical chemists: topological methods are useful for identifying atoms and functional groups on the one hand, and electron pairs, whether bonded or not, on the other. Note that in traditional approaches, the former can only be isolated after projection onto atom-centered basis functions, while the latter can only be made meaningful after arbitrary orbital localization procedures. Broadly speaking, in topological approaches atoms are typically recovered from Bader’s Quantum Theory of Atoms in Molecules (QTAIM), ¹⁶ while electron pairs are recovered from the so-called electron localization function (ELF). ¹⁷ Although the QTAIM is a rigorous theory of quantum mechanical subsystems, while the ELF is one of several heuristic proposals to isolate localization regions from density functions, we will consider them on the same level, since the goal of this introduction is to provide a convenient guide to the use of these two very common tools for the practicing chemist. In this sense, we will give a minimal description of the topological point of view, followed by a basic description of the QTAIM and the ELF as applied to some interesting cases.

Topology

Topological methods in computational chemistry begin, effectively, with the insights of Richard Bader. ¹⁶ These approaches analyze multiple spatial scalar fields, such as the electron density, without external references, providing a robust framework for identifying chemically meaningful regions in real space. Topological analyses partition space exhaustively into disjoint domains and establish well-defined relationships between them. Since the union of the obtained domains reconstructs the whole space, these methods also allow an additive decomposition of quantum mechanical observables, thus bridging the gap between QT and chemical interpretation.

More generally, the concept of a topological space arises from the necessity to generalize the idea of continuity beyond metric spaces, enabling the definition of relationships among abstract objects without relying on a distance function. This fundamental framework is established through a set of elements and a collection of subsets, known as open sets, which must satisfy specific axioms. Very succinctly, a set X together with a collection of subsets U is considered a topological space when they fulfill the following conditions:

1. The empty set ∅ and the entire set X are elements of  U .

2. The intersection of any two members of U is also an element of  U .

3. The union of any number of elements in U remains in  U .

This structure allows spaces to be classified into different topologies. For instance, in a metric space R n , defining neighborhoods using the Euclidean metric naturally leads to the standard metric topology. Additionally, extreme cases include the indiscrete topology, which contains only the empty set and the full set, and the discrete topology, which consists of all possible subsets of X.

To illustrate this, consider a simple set X with two elements: X = a, b. Various possible topologies include:

Each choice defines a unique way to structure the space, affecting its continuity properties and partitioning it into meaningful subsets.

Building upon this foundation, the notion of a neighborhood extends the concept of local surroundings for a given point. A subset N ⊆ X is a neighborhood of x if it contains an open set U such that x ∈ U ⊆ N. This definition generalizes the traditional idea of open intervals in metric spaces.

A function f : X → Y between topological spaces is continuous if the preimage f ⁻¹(U) of every open set U in Y is an open set in X. This aligns with the classical ϵ − δ definition when a metric is introduced. Similarly, closed sets are defined such that a set C is closed if and only if its complement X − C is open. These definitions provide an essential toolkit for analyzing continuity in various mathematical settings.

Beyond static topological considerations, the study of dynamical systems enriches our understanding of how structures evolve. A dynamical system ¹⁸ is characterized by a vector field y defined over an n-dimensional manifold M, governed by differential equations of the form dr/dt = y, where r(t) represents the so-called trajectories of the dynamical system.

In particular, gradient dynamical systems are associated with a potential function ρ : R 3 → R , whose gradient field ∇ρ dictates the system’s behavior:

The integral curves of ∇ρ, also known as flux or field lines, exhibit key properties:

1. Each point in space has a unique trajectory, except at special singular points.

2. The gradient vector is always tangent to the field line at any given point.

3. Field trajectories are orthogonal to isoscalar surfaces, ensuring structured evolution.

4. Trajectories either terminate at a point where ∇ρ = 0 (a critical point) or extend to infinity.

The interplay between topology and dynamical systems offers profound insights into mathematical physics, geometry, and applied sciences. We will summarize some applications to QT in the coming sections, but before that, it is convenient to introduce the concept of critical points.

A critical point (CP) of a scalar function f(x, y, z) is a point (x ₀, y ₀, z ₀) where the gradient of f vanishes:

Critical points are classified into different types based on the behavior of the Hessian matrix at that point:

The eigenvalues of this matrix, λ ₁ ≤ λ ₂ ≤ λ ₃, lead to the following classification:

1. Local Minimum: If λ ₁ > 0, then (x ₀, y ₀, z ₀) is a local minimum.

2. Local Maximum: If λ ₃ < 0, then (x ₀, y ₀, z ₀) is a local minimum.

3. First-order Saddle Point: If λ ₂ > 0 and λ ₁ < 0, then (x ₀, y ₀, z ₀) is a first-order saddle point.

4. Second-order Saddle Point: If λ ₂ < 0 and λ ₃ > 0, then (x ₀, y ₀, z ₀) is a second-order saddle point.

5. Degenerate Case: If any λ _i  = 0, the test is inconclusive.

Critical points play a crucial role in understanding the behavior of scalar fields. Since gradient lines must start and end at CPs (possibly at infinity), each gradient line must connect two of them, known as its α- and ω-limits, respectively. Using a first-order Taylor expansion of the gradient field around a CP, it is easy to show that a local maximum is the ω-limit of a three-dimensional set of gradient lines that together form its attraction basin (see Fig. 1 for a 2D example). Taking into account the previously enumerated properties of gradient lines, it becomes clear that the union of all attraction basins of the maxima of a scalar field fills the space except for a null measure set formed by the 2D (1D) attraction basins of the first and second order saddle CPs, respectively. Thus, a scalar field can be used to partition the entire 3D space into as many regions as there are local maxima: R 3 = ⋃ A Ω A (see Fig. 1-right). The field has induced a topology in R 3 . Attraction basins are separated by local zero-flux surfaces of the gradient field, ∫ _S ∇f ⋅d S  = 0. Another topological partition, the dual topology, can also be induced by taking the α limits of the local minima of the field.

Fig. 1:

The topological partition induced by the f(x, y) = −x ⁴ + 4(x ² − y ²) − 3 field, that displays two local maxima. Left: 3D map of f(x, y). Right: Isolines. The yellow line represents the separatrix (zero flux line) between the two basins in which the system is partitioned.

The electron density and the QTAIM

The electron density is the scalar field that determines the probability density of finding electrons at a given point r , and can be obtained from the system’s wavefunction Ψ as:

where x is an electronic spin-spatial coordinate and we trace out the spin degrees of freedom of the first electron. Within the Born–Oppenheimer approximation, Ψ is parameterized by nuclear positions, allowing the density to be expressed as:

This approximation is sufficient to describe most molecular systems, although thermal effects can be incorporated via statistical mechanics. Most importantly, ρ can be obtained experimentally by X-ray diffraction, Compton scattering, electron diffraction, or neutron scattering. ¹⁹ In X-ray crystallography, for example, accurate diffraction intensities I( k ) allow the reconstruction of ρ( r ),

from the Fourier transform of ρ( r ).

The electron density field exhibits several analytical properties that are relevant in what follows:

1. Nuclear Cusp: According to Kato’s theorem, ²⁰

   (7) − 2 Z α = ∂ ln ρ ¯ ( r ) ∂ | r − R α | r = R α ,

   where ρ ¯ ( r ) is the spherically average density at a given point.

2. Asymptotic Decay: For isolated molecules,

   (8) lim r → ∞ ρ ¯ ( r ) ≃ r 2 ( Z total − N + 1 ) / 2 IP e − 2 r IP ,

   where IP is the first ionization energy.

3. Density Inequalities: The Hoffmann–Ostenhof inequality ²¹ states that

   (9) − 1 2 ∇ 2 ρ + IP − Z r ρ ≤ 0 .

Note that we use atomic units throughout this paper. This means that in Eq. (7), for example, several factors that cancel out the dimensional character of the density are missing. This makes the algebraic derivations easier, but can lead to confusion for the non-expert. Although Kato’s theorem shows that ρ is not differentiable at nuclear positions, a modified homeomorphic field where nuclear cusps are substituted by differentiable maxima allows for a meaningful partitioning of space based on its attraction basins following the topological program sketched in Section 2. This approach underpins the quantum theory of atoms in molecules (QTAIM), that provides an atomic partition (excluding the rare cases where non-nuclear maxima appear in the density) wherein critical points of ρ are associated to chemically meaningful objects. Within the QTAIM, the four types of non-degenerate critical points are named as follows:

1. Nuclear CPs (3, − 3): Local maxima that coincide with atomic positions A. Their attraction basins are called atomic basins Ω _A .

2. Bond CPs (3, − 1): First-order saddle CPs, typically found along the axes that join bonded atoms. Whenever two atomic basins share a common separation surface, a bond critical point (BCP) appears on this surface, and two gradient lines starting at the BCP end at each of the two bonded nuclei. This is called a bond path. The set of bond paths of a system is called the molecular graph, and allows for an automatic determination of the chemical connectivity of a molecule.

3. Ring CPs (3, + 1): Second-order saddle CPs, found normally in cyclic structures at more or less the center of molecular rings.

4. Cage CPs (3, + 3): Minima of ρ that occur within three-dimensional molecular enclosures.

Figure 2 shows the chemical graph of pentaborane(9), B₄H₉. Bond and ring critical points are depicted as small green and red spheres, respectively. The five boron atoms form a square-base pyramid, but interestingly, only the apical boron atom is bonded to the borons of the base, which are linked through hydrogen bridges like those found in diborane. Notice that both the apical to base B-B and the non-bridged B-H bond paths are mostly straight, while the bridging B-H paths are extremely curved inwards. This is characteristic of electron-deficient bonds. In other cases, like in the C-C bonds in cyclopropane, the bond paths curve outwards, signaling annular strain. Each of the four faces of the pyramid that share the apical B atom has a ring critical point, but there is no RCP at the center of the base of the pyramid. Due to this, no cage CP appears in the interior of the pyramid, that we can consider open.

Fig. 2:

Molecular graph of pentaborane(9), computed at the B3LYP//def2-TZVP level. See the text for more details. Graph obtained with the AIMAll suite. ²²

Once the critical points of a system are located, it is possible to compute several other scalar or vector fields at their positions to provide insight into the nature of the chemical bond. For example, the electron density at a BCP, ρ _b , is directly related to the charge accumulation in the bond regions. Since the electron density along a bond axis depends strongly on the nuclear charges of the atoms involved, ρ _b varies significantly when the pair of bonded atoms is changed. Therefore, meaningful comparisons of ρ _b should be limited to systems that share the same bonded atomic species.

Comparative studies indicate that ρ _b correlates well with bond strength. As shown in Table 1, the values of ρ _b for diatomic hydrides provide a clear distinction between bonds traditionally classified as covalent or ionic. Small ρ _b values in LiH and NaH suggest a predominantly ionic character, while larger values in HF and HO indicate stronger covalent interactions. This has led some authors to propose the electron density at BCPs as a quantitative measure of bond order. For example, the ρ _b values for C-C bonds in ethane, ethene, and ethyne are approximately 0.25, 0.36, and 0.43 a.u., respectively, reflecting their classical bond orders of 1, 2, and 3. Extending this approach to benzene yields a calculated bond order n _b  = 1.6, which is consistent with standard chemical intuition. As shown also in Table 1, the relative position of the BCP with respect to the bonded nuclei can be used to define directional bond radii that provide information about the actual size of atoms. Provided that atomic densities decrease exponentially with distance, and that the redistribution of electron density driven by chemical bonding is relatively small in absolute terms, ρ _b is expected to exhibit an exponential dependence on bond distance, reinforcing its role as a fundamental descriptor in chemical bonding theories. The relationship between bond length and bond order, a cornerstone of chemical reasoning, suggests that bond order can be interpreted as a scaled function of ρ _b . Further refinements using second-order density arguments further enhance this understanding.

Several ancillary fields can be constructed from the electron density. One of the most intensively studied is its Laplacian,

Elementary algebra shows that if the sign of ∇² ρ( r ) is positive, the average value of ρ in a neighborhood centered at r is smaller than ρ( r ) itself, and vice versa. This means that the Laplacian of a field monitors its local accumulation (if negative) or depletion (if positive). When computed for an exponential patch of density in a spherical system, the Laplacian shows an alternating sign, with a negative part at small distances followed by a positive one at larger ones. Since in atoms each electronic shell leads to an exponential density patch, the Laplacian provides a way to represent electron shells, and the number of negative and positive regions coincides with the number of shells, at least for atoms with Z ≤ 20. When ∇² ρ is examined in molecules, the inner nuclear shells of atoms are typically not affected by bonding, so they remain as spherical Laplacian regions near the nuclei. The valence shells, however, are strongly distorted by bonding. In many cases the accumulation (∇² ρ < 0) regions of the valence shells of the bonded atoms fuse together, and the Laplacian at the BCP is negative. In other cases, the valence shells replegate toward the nuclei and do not fuse, so ∇² ρ _bcp  > 0. These cases of open and closed shells are characteristic of covalent and non-covalent or ionic interactions, respectively.

The CPs of the Laplacian field have also been extensively studied. Of particular interest are its minima, points of maximum charge concentration (CC). Outside the nuclei and in the valence regions, CCs are found along bond lines, the so-called bonded charge concentrations (BCC), or elsewhere, the non-bonded charge concentrations (NCC), typically associated with lone pairs. Valence CCs tend to appear at distances close to those of the valence shells in free atoms, and it has recently been shown ²³ that they can be associated with either one or both electrons of a Lewis pair.

Figure 3 shows the positions of the BCPs of ρ together with isocontours and the positions of the minima of ∇² ρ in the L i N H 3 + cation. This molecule can be understood as a dative interaction between an ammonia molecule and a Li⁺ cation, similar to ammonia borane. One can easily verify that ∇² ρ has only one nearly spherical electron shell (in red) around the Li nucleus, while it has two around N. The absence of the second shell in Li reinforces its role as a cation, vide infra. There is one Li-N BCP characterized by a positive Laplacian (closed-shell interaction) and three equivalent N-H BCPs with ∇² ρ < 0. It is clear from the isocontours that the valence shells of N and H have fused so that the N-H bonds are of the shared-shell type. Finally, note that there are four CCs around the N atom in a tetrahedral arrangement at about the same distance from the nucleus, so that the spherical valence shell of the free atom has crystallized into four components. The one facing Li is the ammonia lone pair and hosts two electrons. Along each of the N-H directions we also find a CC on the H nucleus. Each one hosts one electron.

Fig. 3:

Some topological descriptors in B3LYP//6-311G* L i N H 3 + . BCPs and the minima of ∇² ρ are shown as small green and yellow spheres, respectively. Isocontours of the Laplacian in a plane containing the Li, N, and one H atom are also shown, with positive(negative) values depicted in blue(red). Notice that the apparent tilt of the plane is due to the perspective. The value of the Laplacian at the BCP, in a.u. is also provided. Graph obtained with the AIMAll suite. ²²

Since R 3 is exhaustively partitioned by these procedures, and quantum mechanical expectation values can be written as integrals over the space, topological approaches provide additive decompositions of expectation values into domain contributions. In the case of one-electron operators, O ̂ , their expectation value becomes a simple sum of domain terms, ⟨ O ̂ ⟩ = ∑ A ⟨ O ̂ ⟩ A . This becomes a double sum of intra- and interdomain contributions in the case of two-electron operators, ⟨ G ̂ ⟩ = ∑ A , B ⟨ G ̂ ⟩ A , B , and so on. ¹² Although the meaning of these domain terms must be carefully examined when the operators involve coordinate derivatives, in many cases the interpretation is straightforward. For example, if we simply use the particle density operator, the number of electrons in a system becomes a sum of domain populations. In the case of QTAIM, each atom has an average fraction of the total electron population, N = ∑ _A N _A . These atomic populations provide orbital invariant atomic charges. In our previous L i N H 3 + example, the QTAIM charges of the Li, N, and H atoms are 0.96, − 1.14, 0.39 a.u., respectively, showcasing that the Li atom is an almost completely ionized species and confirming the expected polarity of the atoms in ammonia from the differences in electronegativity of the N and H atoms.

Using multi-electron particle density operators (see next section for pair density), one can also partition the total number of electron pairs, trios or, more generally, electron n tuples into two, three or n regions. If we consider the electron population of a region as a random variable, we can obtain e.g. its variance in a region or the covariance between two regions from the pair partition. In recent years it has become clear that bonding is due to mutual correlation between electron populations in atomic domains. Thus, if a change in the population of one atom simultaneously affects the population of n − 1 other atomic domains, we can speak of an n-center bond. In the two-center case, the covariance of the electron populations is a direct measure of their bonding order. It is usually transformed into the delocalization index, DI(A, B) = −2cov(n _A , n _B ). ²⁴ In L i N H 3 + , the calculated DIs are 0.104 and 0.852 for the Li-N and N-H interactions. The former is a very weak interaction, while the latter is a polar covalent single bond. DIs exist between atoms that are not directly bonded and can be used to measure the strength of non-covalent interactions. They are 0.002, 0.015 for the Li-N and H-H interactions.

These ideas can be taken several steps further, and the full probability distribution of finding an integer number of electrons in a set of domains can be defined and computed. These are called electron distribution functions (EDFs). ²⁵ As an example, consider a set of three atoms (or groups of atoms) in a molecule. Then we can obtain p(n ₁, n ₂, n ₃), the probability that n ₁ electrons are found in the region Ω₁, n ₂ in Ω₂, n ₃ in Ω₃, and the N − (n ₁ + n ₂ + n ₃) remaining electrons in R 3 − Ω 1 ∪ Ω 2 ∪ Ω 3 , the complementary domain, for all sets of n ₁, n ₂, n ₃ integers such that n ₁ + n ₂ + n ₃ ≤ N. It can be shown that from these probabilities, which are close to Pauling resonance structures, all multicenter bond orders can be obtained. Of course, we can divide each of these n _i populations into up- and down-spin subsets and ask for the probability that n 1 α and n 1 β electrons are in Ω₁.

Coming back to L i N H 3 + , we show two differently grained EDF descriptions. First we divide the system into the Li and NH₃ moieties. The most likely arrangement of the 12 electrons of the system is that with two electrons in the Li domain and ten in the NH₃ one: p(2, 10) = 0.946. This is followed by p(3, 9) = 0.047. A very minor p(1, 11) = 0.007 value is also found. This distribution shows that the lithium is a cation, and that an electron from the lone pair of the nitrogen is involved in an exchange interaction with it less than 1 % of the time. A second possibility is to explore the full distribution, considering the five atoms of the system at the same time. For the order Li, N, H1, H2, H3 the largest probability is p(2, 10, 0, 0, 0) = 0.109. This reflects the highly polar nature of the N-H bond. The next term is p(2, 9, 1, 0, 0) = 0.094, where one of the H atoms holds one electron and that is obviously triply degenerate. Next come p(2, 8, 1, 1, 0) = 0.080, p(2, 7, 1, 1, 1) = 0.065, etc. We will stop the discussion here, but a full analysis of a system’s EDF provides a great deal of relevant information.

The electron localization function

The Electron Localization Function (ELF) was originally introduced by Becke and Edgecombe ¹⁷ to describe the likelihood of finding an electron near another reference electron of the same spin. It is derived from the conditional probability of finding an electron at a point r ₁ when another is known to be located at r ₂. The two-electron probability distribution function is expressed as:

Recall that x is a spin-spatial electron coordinate. This function captures electron correlation effects, reflecting both the Pauli exclusion principle and Coulomb repulsion. In the absence of correlation, the probability simplifies to a product of independent electron densities:

where normalization ensures that if one electron is known to be at x ₁ there are only N − 1 remaining electrons to be located at x ₂. If we consider two electrons with the same spin, e.g. up-spin, then ρ 2 α α ( r 1 , r 2 ) must vanish if r ₁ =  r ₂ due to Pauli exclusion. Then, if r ₂ is close to r ₁, r ₂ =  r ₁ +  s , a Taylor expansion with respect to s = | s | of the spherically averaged value of this same-spin probability must start with a quadratic term 1 2 a s 2 , where a is the so-called curvature of the Fermi hole. If a is small, then there is a large region around a reference electron located at r ₁ where no second same-spin electron can be located, and the reference electron is said to be localized. The ELF measures the localization of electrons by analyzing the curvature of the Fermi hole, which can be obtained at the single-determinant level:

Since D ^α depends on the density, it is not an intensive quantity and its values are difficult to compare across different systems. To standardize ELF values, normalization to the homogeneous electron gas is introduced by computing D ^α for an electron gas with the same density as the system under study:

The ELF kernel is then defined as:

and a final Lorentzian transformation maps ELF to the fixed range [0, 1]:

ELF values close to 1 indicate strong electron localization, while values near 0.5 correspond to the homogeneous electron gas behavior. The original ELF definition was restricted to Hartree–Fock (HF) wavefunctions, restricting considerably its scope. Savin et al. ²⁶ extended the ELF to density functional theory (DFT), utilizing the positive definite kinetic energy density, t( r ) = (1/2)∇ ⋅∇′ρ( r ′; r ), where ρ( r ′; r ) is the first order reduced density matrix. It can be shown that:

where the second term is the kinetic energy density of a system of bosons with the same density as that of the electronic system. ²⁷ From this, the Pauli kinetic energy density, representing fermionic contributions, is defined as:

Normalization using the Thomas-Fermi kinetic energy density gives:

which are the equivalents to Eqs. (15) and (16). This formulation provides a more general framework for ELF calculations within DFT, ensuring a broader applicability to big systems and condensed matter.

As a scalar function, the ELF can be subjected to topological analysis, its critical points can be obtained, and R 3 can be partitioned into an exhaustive set of non-overlapping attraction basins associated with its maxima. It is shown that these localization regions correspond exactly to those associated in chemistry with electron pairs: nuclei, bonds, lone pairs. Other more exotic cases are also found, mapped for instance to 3-center-2 electron bonds, localized isolated electrons in electrides, etc. Of course, any type of expectation can also be decomposed into basin contributions, and ELF basin populations, variances, or covariances can be reported.

In order to illustrate the utility of ELF, we start with atoms. ELF effectively recovers atomic shell structures in agreement with the Aufbau principle, ²⁸ with as many maxima as shells. The elements of the first row are shown in Table 2 to illustrate it. If we integrate the electron populations over these domains we check that they correspond to the expected 2n ² behavior. It has also been shown that core populations remain nearly constant across different molecular environments, supporting the concept that core electrons are largely unaffected by chemical bonding.

Moving from atoms to molecules, ELF provides insights into electron localization in chemical bonds and lone pairs. ²⁹ The synaptic order of an ELF basin, defined as the number of basins directly connected (sharing a first-order saddle) to it, classifies electron localization as follows: ³⁰

1. Monosynaptic basins (M) correspond to lone pairs and often exceed two electrons (typically due to higher electronegativies of the atom holding the lone pair).

2. Disynaptic basins (D) correspond to two-center covalent bonds and typically hold slightly less than two electrons. For instance, in the CH₃X series (X = C, N, O, F), bond populations decrease as electronegativity increases: N(Ω) = 1.82 (C–C), 1.57 (C–N), 1.22 (C–O), 0.85 (C–F), reflecting increased polarization. Note that Protonated disynaptic basins show electron populations between 1.3 and 2.5 a.u., consistent with the VSEPR theory.

3. Trisynaptic or polysynaptic basins (T) describe multi-center bonding.

For example, in water, the ELF analysis reveals two distinct monosynaptic maxima for the lone pairs forming a distorted tetrahedron (Fig. 4a) whereas for N₂, lone pairs appear at the external extremes of the N atoms (Fig. 4b). Moreover, the ELF can be used to distinguish bond multiplicities based on the number and shape of disynaptic bonding attractors. Ethane (C₂H₆) has a single ELF attractor (Fig. 4c) for the C–C bond, whereas ethylene (C₂H₄) exhibits two (Fig. 4d).

Fig. 4:

ELF isosurface = 0.9 for (a) water, (b) N₂, (c) ethane, (d) ethylene.

Note that the topology of the ELF also validates the VSEPR model ³¹ by providing quantum mechanical confirmation of the spatial distribution of electron pairs. For molecules of type AX ₄, ELF localization domains align with predicted electron pair geometries, such as:

1. Tetrahedral arrangements in AX ₄ molecules.

2. Trigonal bipyramidal lone pair distributions in AX ₅ molecules.

3. Cylindrical lone pair structures in AX ₆ molecules.

ELF basins are also useful for analyzing electronic fluctuations. The population variance σ ²(Ω) within a given basin is calculated as:

where ρ ₂(r ₁, r ₂) is the spinless pair function obtained after integrating the spin coordinates in Eq. (12). To facilitate comparisons, a relative fluctuation(λ(Ω)) can be defined as:

The study of ELF fluctuations in conjugated organic systems has provided interesting insights into the role of resonance and the effect of electron correlation in molecules. For example, a study in oxaziridine has shown small and fluctuating electron densities in the ELF valence basin between the nitrogen and oxygen atoms (N-O bond) that highlights the correlated nature of these electrons. ³²

Finally, a point on method dependency. Table 3 shows the effect of increasing the basis set size on ELF populations for H₂O. Variations are minor when different basis sets are used, hence, this is not a major point of concern when applying in general topological methods.

Combined QTAIM-ELF-EDF analysis

The study of molecular structure requires a comprehensive approach that integrates multiple theoretical frameworks. The Quantum Theory of Atoms in Molecules (QTAIM), the Electron Localization Function (ELF), and other available real space procedures provide complementary insights into bonding that may be leveraged together. Here we apply this combined approach to a very simple, yet intriguing system: the carbon suboxide molecule, C₃O₂. This is a well known stable oxide of carbon. ³³ Initial spectroscopic evidence supported a linear structure in the gas phase, later corroborated by simple calculations, leading to the cumulenic structure O=C=C=C=O usually found in textbooks. However, better resolution infrared experiments, ³⁴ showed that the molecule had a very labile bent structure, with C–C–C angle of 156°. This is corroborated by high level coupled-cluster results, and has led Frenking and coworkers to introduce the so-called carbones, ³⁵ where a C(0) central moiety binds by dative bonds to CO molecules, leaving lone pairs that bend the system, OC → C ← CO. Although standard density functionals fail to provide a bent structure, a simple Hartree–Fock calculation leads to an angle close to 140°, that opens to 159°, very close to the experimental value, if an MP2/def2-TZV optimization is performed. To avoid problems with the meaning of bonding descriptors like the delocalization index or the ELF in the case of perturbative calculations, while allowing for correlation effects to be reasonably taken into account, we will just use the MP2 geometry and report M06-2X/def2-TZV DFT results to illustrate the use of combined analyses.

We first note that, at the same level of theory, the CO molecule has QTAIM atoms with charges equal to ±0.929 a.u., a density at the BCP ρ = 0.421 a.u., and a DI equal to 2.02 (polar multiple bonds typically display DIs smaller than the nominal values, a triple bond in this case). The Laplacian at the BCP is only slightly negative, ∇² ρ = −0.017 a.u., but the energy density H = −0.692 a.u. is clearly negative. The ELF V(C, O) basin hosts 2.34 electrons, and the V(C) and V(O) lone pairs, 2.74 and 4.72 e, respectively.

QTAIM analysis

Figure 5 shows that upon formation of the C₃O₂ suboxide, the valence shell of the central C1 is fused with those of its C2 and C3 neighbors. Actually, the value of ∇² ρ at the C1-C2 BCP is −0.792 a.u., and H = −0.394 a.u., with a C1–C2 DI of 1.712. Upon complexation with the C atom, the CO molecules suffer several noticeable changes in their electronic structure. Their bond length increases from 1.168 to 1.204 Å, and ρ and ∇² ρ at the BCP change to 0.397, and −0.597 a.u., respectively, with H = −0.651. Similarly, the DI decreases to 1.71. In the suboxide, the charge of each CO molecule is −0.030 a.u., and Q(O) = −0.801 a.u., leading to Q(C1) = 0.059 a.u. All of these findings are compatible with the expected bonding mode of a classical carbonyl ligand with π-backbonding interactions that weaken the CO bond. Notice that the DIs support a cumulene-like picture with strong C₁=C₂ links.

Fig. 5:

The carbon suboxide system. Top: Bond critical points of ρ (green spheres), valence minima of ∇² ρ (yellow spheres), bond paths (black lines), and projection of the interatomic surfaces on the molecular plane. Isocontours of the Laplacian are also shown in red and blue for negative and positive values, respectively. Bottom: η = 0.8 isosurfaces of the ELF, color-coded according to attractor types. Disynaptic in green, monosynaptic in red.

On top of this description, the position of the minima of the Laplacian (yellow spheres in the top panel of Fig. 5) shows that a notorious non-bonded charge concentration exists along the C ₂ axis that would support the carbone model. Interestingly, each of the C1–C2, C1–C3 bonds shows a pair of CCs along the bond path. The ones belonging to the C1 basin are extremely close to the C1–C2 and C1–C3 interatomic surface (pink lines). This is also in agreement with a donor-acceptor interaction, where the two-electrons of the σ components of these links would come from the donating carbonyls. Interestingly, the three-center delocalization index (a measure of the intensity of 3-center interactions) between the central carbon and the two carbonyls is very large, equal to 0.393.

EDF analysis

A convenient partition of the carbon suboxide into fragments is that of two carbonyls and the central carbon atom. In this three-domain partition, which we will order as (CO1, C, CO2), the largest probability electron distribution is p(14, 6, 14) = 0.128, followed by p(13, 7, 14) = p(14, 7, 13) = 0.087, p(14, 5, 15) = p(15, 5, 14) = 0.096, and p(13, 6, 16) = p(15, 6, 13) = 0.067. This supports a small positively charged central carbon engaging in covalent exchanges with the carbonyls.

If we isolate the C1 atom, the probability p(n _C1) is found in the top panel of Fig. 6. The computed values, are non-negligible for electron counts ranging from n = 2 to n = 10. This indicates that the C1 atom is involved in four different links, with a final distribution slanted toward n < 6 counts that reveal a slightly positive charge, as already put forward. This is corroborated if we also include the theoretical binomial distribution of four independent symmetric bonds, p ( 2 + n ) = 8 n / 2 8 , that is plotted as a continuous green line. If we average the n _C1 = 6 ± m probabilities, where m runs from 1 to 4, as found in the filled cyan circles, the match between the modeled and computed results is close to perfect. This shows that even the lone pair is participating in bonding with the carbonyls, and that the electronic structure of the suboxide contains elements of both the cumulene and carbone descriptions. This picture can be complemented by the bottom panel, which shows that each of the carbonlys engage in two links with the rest of the system. The match between the average probabilities and the binomial model for a double bond is, as in the C1 case, very good. Notice that the distribution is slightly skewed towards greater than nominal electron populations.

Fig. 6:

Top: Probability of finding a given number of electrons in the C1 basin. The open triangles are the computed values, the green line is the binomial distribution corresponding to four independent symmetric covalent bonds, and the cyan full circles the average of p(6 + n) and p(6 − n) for n = 1, …, 4. Bottom: Equivalent plot for one of the carbonyl. Now the number of bonds of the model is two.

This worked example has demonstrated how the Electron Localization Function (ELF) complements the Quantum Theory of Atoms in Molecules (QTAIM) by providing insights into bonding and lone pair interactions in the carbon suboxide (C₃O₂) moiety, an angular molecule with a low bending force constant. The central carbon’s lone pairs have been linked to the concept of carbones, a family of zero-valent compounds extended to other group 14 elements. While the Laplacian of the electron density suggests lone pair and donor-acceptor interactions, the ELF analysis supports the carbone model. However, delocalization indices (DI) and electron density functions (EDFs) indicate four independent bonding channels, reinforcing a traditional cumulenic structure. These findings highlight the importance of combining ELF and QTAIM for a comprehensive understanding of molecular bonding, especially in unconventional bonding scenarios.

Summary and conclusions

Quantum mechanics has profoundly influenced chemistry, yet fundamental conceptual differences persist between the methods and language used in both disciplines. While quantum theory (QT) provides a comprehensive theoretical foundation, chemical intuition remains rooted in empirical models. Efforts to integrate QT with chemical reasoning, particularly through topological analysis, have yielded promising results.

We have briefly reviewed two main methodologies within this approach, the quantum theory of atoms in molecules (QTAIM) and the electron localization function (ELF). The topological analysis of electron density proposed by the QTAIM offers a rigorous, quantitative tool for defining atoms within a molecule. Notably, the electron density at bond critical points (BCPs), denoted as ρ _b , serves as a powerful indicator for characterizing bond strength and bond order. Once an atomic partition of the space has been obtained, all sort of observables can be additively decomposed. We have revised delocalization measures as the delocalization index (DI), and the probabilities of finding a given distribution of the electrons of the system in the atomic regions, the so-called electron distribution functions (EDFs). On the other hand, the Electron Localization Function (ELF) provides an intuitive and quantitative measure of electron localization and electron pairing effects, enabling the identification of covalent bonds and lone pairs that match the Valence Shell Electron Pair Repulsion (VSEPR) theory.

The combined use of the ELF and the QTAIM offers a comprehensive framework for analyzing the electron distribution in molecular systems, providing important insights into charge distribution, bond character, and electron correlation effects. As an example, we showed in this contribution how ELF can complement QTAIM by revealing details about bonding and lone pair interactions in the carbon suboxide molecule, C₃O₂, which is an angular molecule with a very low bending force constant. The angular geometry suggests the existence of lone pairs at the central C atom, and this has been used to propose a family of zero-valent compounds called carbones that have been generalized to other elements of group 14 such as germanium. As shown, a single approach may give rise to partial, even misleading conclusions. For instance, the position of the charge concentrations in the Laplacian of the density suggests a lone pair and dative donor-acceptor interactions. The ELF image also clearly suggests that the carbone picture should be favored. However, the DI values and the EDFs also indicate that the C1 engages in four independent delocalization channels (bonds), in agreement with the traditional cumulenic structure for this molecule. These findings emphasize the significance of a combined approach, utilizing both ELF and QTAIM, for a deeper understanding of the chemical bonds in complex molecular systems. Further exploration of these methods, particularly in systems with unusual bonding scenarios are welcome. We firmly believe that topological analyses provide a robust approach to rationalize the electronic structure of systems and chemical processes. Unlike other methods, they remain largely independent of the computational basis set and are easily correlated with chemical intuition.