Chemical bonding concepts emerge naturally from maximally entangled atomic orbitals

The chemical bond is a central concept in chemistry, originating from Gilbert N. Lewis’s idea of electron pairs shared between atoms¹. Lewis’s theory, developed independently of quantum mechanics, is remarkably simple and profoundly influential. It is thus commonly taught already at the high school level. Heitler and London provided the first quantum mechanical description of a hydrogen molecule bond², later evolving into valence bond theory^3,4, which uses superpositions of configurations of electrons occupying atomic orbitals. Molecular orbital theory⁵ followed historically, combining atomic into molecular orbitals before occupying them with electrons. Both approaches link bonding to quantum mechanics by virtue of rudimentary wave functions, guided by chemical intuition, and together they form the foundation of modern bonding theories.

Contrarily, modern electronic structure methods based on sophisticated wave function ansätze offer great accuracy^6,7,8,9,10 but at the same time complicate intuitive bonding descriptions. To address this issue, various chemical bonding tools have emerged^{11,12,13,14,15}, categorized into real-space or topological, and Hilbert space approaches. Concerning the former category, approaches such as Bader’s quantum theory of atoms in molecules¹⁶ and the electron localization function^17,18 rely on topological analyses of the electron or same-spin pair density. Consequently, they require computationally intensive integrations over three-dimensional regions to derive descriptors like bond order and aromaticity indices^19,20. While being largely basis-set independent, these methods therefore suffer from scalability issues and integration errors²¹. Due to these issues, a recent information theoretical approach to bonding analysis based on real-space partitioning, though conceptually interesting, has thus far been limited to mean-field descriptions²². In contrast, Hilbert space methods^11,23,24,25 partition molecular systems into atomic domains using basis functions, offering computational efficiency but facing challenges in basis set dependence and ambiguity in the orbital assignments²⁶. Recent orbital localization schemes, such as intrinsic atomic orbitals (IAOs)²⁵, mitigate these limitations, providing a robust foundation for chemical bonding analyses, yet without inherently yielding bonding concepts.

Despite these challenges, Hilbert space partitioning methods hold significant promise. As we will demonstrate, these methods are particularly well-suited for applying tools from quantum information theory (QIT) to analyze chemical bonding, offering a perspective that addresses their current limitations. To further motivate our work, we recall that recent studies^{27,28,29,30,31,32,33,34,35,36,37} have applied QIT tools to partitions based on molecular orbitals, thus quantifying electron correlation and the representational complexity of molecular quantum states. In contrast, quantifying bonding structures with an emphasis on bonding orders is a significantly more challenging task, as it requires adapting the QIT approach to the nonlocal nature of covalent bonds by instead referring to fully localized orbitals. The primary accomplishment of our work is to devise a general scheme for computing maximally entangled atomic orbitals (MEAOs), whose entanglement patterns quantitatively recover both Lewis (two-center) and beyond-Lewis (multicenter) bonding structures, with multipartite entanglement serving as a robust index of bond strength.

Accordingly, by leveraging parallels between chemical bonding and quantum entanglement, our work provides a QIT-based framework that captures crucial features of chemical bonding such as hybridization, bond orders, multicenter bonding, conjugation, and aromaticity without a priori chemical assumptions. Our framework, which we validate on standard and unconventional bonding cases, thus provides a new framework for understanding fuzzy chemical concepts using rigorous, quantitative descriptors from quantum information theory.

Entanglement in a single covalent bond

We motivate our general QIT-based framework by first analyzing the idealized covalent bond formed by two electrons in a symmetric diatomic molecule. For this, let φ_L and φ_R be the two relevant (real-valued) atomic orbitals localized on the left and right atomic centers, with an overlap \(S=\left\langle {\varphi }_{{{{\rm{L}}}}}| {\varphi }_{{{{\rm{R}}}}}\right\rangle \ne 0\). The prototypical bonding state \(\left|{\Psi }_{{{{\rm{bond}}}}}\right\rangle\) is then obtained by occupying the corresponding bonding orbital with the electron pair,

where \({\phi }\) and \(\bar{\phi }\) denote the bonding and antibonding orbitals,

Here and in the following, we use the formalism of ‘second quantization’ and recall that any orbital φ gives rise to a four-dimensional Fock space, spanned by the states \({\left|0\right\rangle }_{\varphi },{\left|\uparrow \right\rangle }_{\varphi },{\left|\downarrow \right\rangle }_{\varphi }\), and \({\left|\uparrow \downarrow \right\rangle }_{\varphi }\), where the arrow denotes the electron spin and \({\left|0\right\rangle }_{\varphi }\) the vacuum state. It is crucial to note that the idealized state, Eq. (1), takes the form of a product state (in 2^nd quantization) and, accordingly, it does not contain any entanglement between the bonding and antibonding orbitals. This is by no means a contradiction to the common expectation: Entanglement is a relative concept and depends on the division of the total system into orbital subsystems. By referring to molecular orbitals \(\phi,\bar{\phi }\), the corresponding orbital entanglement merely quantifies the validity of the independent electron-pair approximation rather than the bonding order. This, in turn, reflects very well the conceptually different perspectives of molecular orbital and valence bond theory.

According to the definition of effective bond order (EBO)¹⁵, which is the difference between the occupation numbers of the bonding and antibonding orbitals divided by 2, the state in Eq. (1) represents a “perfect” single bond with EBO = 1. In general, to compute the EBO for more realistic molecular wave functions, one would first need to categorize various molecular orbitals as bonding, antibonding, or nonbonding. This can be based on factors such as spatial symmetry, or their ability to promote or inhibit electron sharing between the atomic centers^11,38. However, such criteria can become ambiguous and arbitrary in multicenter molecules, potentially leading to erroneous results. As our work will show, it is precisely the QIT framework used in a valence bond theoretical context that offers excellent prospects for overcoming these deficiencies of approaches based on molecular orbital theory.

To elaborate further on the idealized covalent bond and align our QIT perspective with its nonlocal character, we introduce the symmetrically orthogonalized atomic orbitals

In this orbital basis, the state in Eq. (1) has the following form:

That is, it is maximally entangled relative to the orbitals \({\widetilde{\varphi }}_{{{{\rm{L}}}}}\) and \({\widetilde{\varphi }}_{{{{\rm{R}}}}}\), with an entanglement value \(E=\log (4)\). Here, we used the fact that the entanglement E for pure states follows as the von Neumann entropy \(S(\widehat{\rho })\equiv -{{{\rm{Tr}}}}[\widehat{\rho }\log (\widehat{\rho })]\)³⁹,

of the corresponding orbital reduced density matrices (RDM)

We also recall that the maximal entanglement between two subsystems of dimension d (in our case d = 4) is simply \({E}_{\max }=\log (d)\). This entanglement value indicates maximal correlation between the physical observables measured on the left and right orbital, which, by construction, recovers crucial bonding features such as electron sharing and spin pairing.

Finally, we would like to stress that the orbitals \({\widetilde{\varphi }}_{{{{\rm{L}}}}}\) and \({\widetilde{\varphi }}_{{{{\rm{R}}}}}\) define a partition of the underlying one-particle Hilbert space, which conceptually resembles the real-space partitioning in the quantum theory of atoms in molecules¹⁶. In more general situations involving multiple atomic centers, or when each atomic center hosts more than one atomic orbital, calculating the entanglement between any two atomic subspaces becomes significantly more challenging for several reasons. First, the reduced state \({\widehat{\rho }}_{AB}\), defined on two atomic centers A and B, is high-dimensional because the Fock spaces of A and B grow exponentially with the number of orbitals on each center. Second, \({\widehat{\rho }}_{AB}\) is generally a mixed state as a result of the coupling of A and B to other atomic centers. Consequently, the closed formula Eq. (5) for entanglement in pure states is no longer applicable, and the entanglement between A and B must instead be determined numerically as the minimum “distance” from \({\widehat{\rho }}_{AB}\) to the set of unentangled mixed states S_sep^40,41

Nevertheless, the high dimensionality of \({\widehat{\rho }}_{AB}\) can often be circumvented when the entanglement between A and B is predominantly captured by a few pairs of localized orbitals \({\widehat{\rho }}_{ij}\). This is a crucial aspect of our work which we will elaborate on in the next section.

Maximally entangled atomic orbitals

The centerpiece of our results is the set of maximally entangled atomic orbitals (MEAOs). We briefly outline their construction and usage here, while detailed procedures are provided in the Methods section. Given a Hilbert space partition \(\Pi=\{{{{{\mathcal{H}}}}}^{(1)},{{{{\mathcal{H}}}}}^{(2)},\ldots,{H}^{(M)}\}\) with atomic subspaces \({{{{\mathcal{H}}}}}^{(m)}\), MEAOs are a basis of localized orbitals \({{{{\mathcal{B}}}}}^{{{{\rm{MEAO}}}}}=\{\{{\phi }_{i}^{(1)}\},\{{\phi }_{j}^{(2)}\},\ldots,\{{\phi }_{k}^{(M)}\}\}\) that ideally maximizes the total inter-center entanglement, defined as the sum of entanglement between two orbitals sitting on different atoms

where \({\widehat{\rho }}_{ij}\) is the orbital RDM of orbital i and orbital j that sits on different atomic centers, obtained by tracing out all other orbital degrees of freedom from the overall wave function \(\left|\Psi \right\rangle\). While Eq. (8) is physically well motivated, its direct evaluation is computationally demanding due to the mixed-state nature of \({\widehat{\rho }}_{ij}\) and its dependence on the four-particle RDM. To identify the MEAOs, we need a computationally affordable objective function to guide the orbital rotations. To this end, we observe that the entanglement in a two-orbital RDM \({\widehat{\rho }}_{ij}\) is maximized when it exactly matches the single covalent bond state, Eq. (4). The entanglement in the state, Eq. (4), can be attributed to two important superpositions: (i) \({\left|0\right\rangle }_{{{i}}}\otimes {\left|\uparrow \downarrow \right\rangle }_{{{j}}}\) and \({\left|\uparrow \downarrow \right\rangle }_{{{i}}}\otimes {\left|0\right\rangle }_{{{j}}}\); (ii) \({\left|\uparrow \right\rangle }_{{{i}}}\otimes {\left|\downarrow \right\rangle }_{{{j}}}\) and \({\left|\downarrow \right\rangle }_{{{i}}}\otimes {\left|\uparrow \right\rangle }_{{{j}}}\). These two superpositions are purely two-body, and the coherence therein is measured by the elements \({\Gamma }_{j\uparrow,j\downarrow }^{i\uparrow,i\downarrow }\) and \({\Gamma }_{j\uparrow,i\downarrow }^{i\uparrow,j\downarrow }\), respectively, of the two-particle reduced density matrix (2RDM) \({\Gamma }_{k\sigma,l\tau }^{i\sigma,j\tau }=\left\langle \Psi \right|{f}_{i\sigma }^{{{\dagger}} }{f}_{j\tau }^{{{\dagger}} }{f}_{l\tau }{f}_{k\sigma }\left|\Psi \right\rangle\). Based on this observation, we define the objective function for obtaining the MEAOs as

where orbitals i and j belong to different atomic centers, and

where the minimization is performed over all bases of localized orbitals that respect the Hilbert space partition Π (denoted as \({{{\mathcal{B}}}} \sim \Pi\)). In Fig. 1, we summarize the workflow of using MEAO for bonding pattern analysis:

1. a.

   Choose a Hilbert space partitioning.

2. b.

   Construct MEAOs using a mean field wave function \(\left|{\Psi }^{{{{\rm{MF}}}}}\right\rangle\) and Eqs. (14), (10).

3. c.

   Compute the mutual information \({I}_{ij}=S({\widehat{\rho }}_{i})+S({\widehat{\rho }}_{j})-S({\widehat{\rho }}_{ij})\) for all pairs of inter-center orbitals i and j, and apply graph-theoretical methods to group the MEAOs into bonding pairs or bonding groups.

4. d.

   Evaluate the pairwise or multipartite entanglement (defined in section “Beyond-Lewis structures: genuine multipartite entanglement”) for these MEAO groups, using a correlated wave function \(\left|\Psi \right\rangle\).

The alternative text for this image may have been generated using AI.

Full size image

a Perform a Hilbert space partition with localized orbitals. Each group of localized orbitals defines the local Hilbert space for each atom. b Construct MEAOs via internal orbital rotations that preserve the local Hilbert spaces, and minimize the objective function, Eq. (14), computed with a mean-field wave function. c Group MEAOs into two-center and multicenter bonds by identifying strong links in the orbital correlation graph. MEAOs are represented by the vertices in the correlation graph whose colors indicates the atomic assignment, whereas strong and weak orbital correlations are illustrated with the solid and dashed edges, respectively. d Apply quantum information theoretical tools to each identified bond.

Lewis structures and bipartite entanglement

Recovering standard chemical concepts: Within the framework of standard Lewis theory, molecular bonding is described by Lewis structures, where bonds are represented by pairs of shared electrons depicted as lines between atoms. The bond multiplicity is then given by the number of lines connecting two atoms. What sets MEAOs apart from other localized orbitals is their ability to naturally reflect the Lewis structures of molecules through their correlation and entanglement properties.

In Fig. 2, we compare the shapes and orbital correlations of the valence IAOs and MEAOs of ethene (C₂H₄) with the cc-pVDZ basis. The MEAOs are constructed from the 14 minimal IAOs using a Hartree-Fock solution. Subsequently, two complete active space (CAS) calculations are performed using the IAOs and MEAOs, respectively, to compute orbital correlation and entanglement. While the IAOs retain much of their free atomic character, with a clear distinction between one 2s and three 2p orbitals, the MEAOs exhibit the classical sp² hybridization of carbon, forming three sp² orbitals and one 2p orbital. This remarkable emergence of hybridization—a fundamental concept in organic chemistry—demonstrates how MEAOs naturally capture chemical intuition (see Supplementary Note 2 for further examples).

The alternative text for this image may have been generated using AI.

Full size image

a Isosurfaces of valence carbon IAOs. b Isosurfaces of valence carbon MEAOs. c Orbital correlation graph of IAOs. Values of the normalized single orbital entropy \(S({\widehat{\rho }}_{i})/\log (4)\) and the normalized orbital-orbital correlation \({I}_{ij}/\log (16)\) are represented in log-scale by the color of the nodes and edges of the graph, respectively. d Same as c but for MEAOs. The calculations are performed with the cc-pVDZ basis.

The differences in orbital shapes between IAOs and MEAOs result in distinct correlation graphs, where MEAOs provide a sparse and interpretable picture of bonding. High pairwise correlation values (\({I}_{ij}\approx \,{I}_{\max }\)) clearly identify bonding interactions, with non-bonding pairs showing negligible correlation values (\({I}_{ij} \sim 1{0}^{-3}{I}_{\max }\)). For ethene, the correlation graph of the MEAOs reveals six bonding interactions: four C–H σ bonds and two C-C bonds (one σ and one π), as shown in Table 1. In contrast, the correlation graph of the IAOs is dense and does not provide a clear bonding picture.

For each two-center bond, there exists a pair of MEAOs with correlation (\({I}_{ij}/{I}_{\max }\)) and entanglement (\({E}_{ij}/{E}_{\max }\)) values close to 1. The number of such pairs matches the bond orders of the molecules analyzed in Table 1. Furthermore, deviations of \({I}_{ij}/{I}_{\max }\) or \({E}_{ij}/{E}_{\max }\) from 1 reflect deviations of the ground state from the idealized Lewis structure⁴². The results also highlight chemical trends, such as that π bonds being generally lower in covalency than σ-bonds. In summary, MEAOs reorganize the ground state wave function into a representation as closely aligned as possible with the Lewis structure of the molecule, providing a robust framework for quantifying bond orders through orbital-orbital entanglement.

Before we move on to more advanced applications of the MEAO formalism, we first point out that orbital entanglement directly captures only bonding effects arising from electron sharing and spin coupling, namely covalent bonding. Ionic or intermolecular bonding effects are only indirectly shown in the reduction of orbital entanglement. To demonstrate this constraint, we studied three molecules with increasing levels of ionicity, Li₂, LiH, and LiF. As shown in the second panel in Table 1, we indeed observe a consistent decrease going from Li₂ to LiF, in both mutual information and entanglement between the two most correlated MEAOs. Finally, we test against the van der Waals molecule He₂, which was found to have zero correlation and entanglement between the MEAOs on the two helium atoms. We remark that this number of absolute zero is ensured by the fact that the minimal IAO subspace from which the He-MEAOs are constructed consists of only one orbital for each atom. And since both electrons are placed on this orbital, the total ground state is simply a product state between the two fully occupied He-MEAOs. To fully confirm the lack of covalency in He₂, we extended the minimal IAO subspace to include the full basis set, only to find a vanishingly small amount of correlation between the two atomic centers.

The challenging harpoon mechanism: While the singlet ground state of LiH at equilibrium (R_Li-H ≈ 1.6 Å) is predominantly ionic, it transitions to a covalent bonding character of the singlet first excited state at an avoided crossing around R_Li-H = 3 Å (see Fig. 3a). This transition is accompanied by a shift of electron density from the hydrogen atom towards the center of the molecule. Therefore, when the molecule is stretched to dissociation, the electron sharing first increases as it enters the covalent phase, and then decreases due to final dissociation. This process is described as the harpoon mechanism⁴³. The hallmark of this mechanism is a peak in the covalent bond order around the avoided crossing⁴⁴, which is confirmed by the electron delocalization index δ_Li-H^20,42 (Fig. 3b). This index measures the covariance of electron populations in the Li and H atoms within a real-space partition (see Supplementary Note 4 for a definition). The bonding in this system is particularly challenging to describe because standard Mulliken-like electron-sharing analyses based on Hilbert space partitioning fail to detect the signature peak in the bond order⁴⁵. By contrast, the quantum information approach using MEAOs successfully identifies this challenging feature, underscoring the strength of our framework in capturing complex bonding phenomena.

The alternative text for this image may have been generated using AI.

Full size image

a Energies (in atomic units, a.u.) of the lowest two states in the spin singlet sector (X¹Σ⁺ and A¹Σ⁺) and the triply degenerate first excited state in spin triplet sector (X³Σ⁺). b The electron delocalization index δ_Li-H of the singlet ground state. c The highest entanglement value (normalized by \(\log (4)\)) between two maximally entangled atomic orbitals, one localized on Li and and one on H, in the thermal state at β = 10³ Ha⁻¹ involving the singlet ground state and the triply degenerate triplet first excited states.

In Fig. 3, we present the low-lying energy spectrum \(\{{{{{\mathcal{E}}}}}_{i}\}\) of LiH calculated using the aug-cc-pVDZ basis set, the electron delocalization index δ_Li-H of the ground state (taken from ref. ⁴⁴), and the highest entanglement between a Li-MEAO and a H-MEAO in the approximate thermal state. This thermal state includes the lowest four eigenstates \(\left|{\Psi }_{i}\right\rangle\), given by

where β = 10³ Ha⁻¹ represents the inverse temperature. Incorporating a low-temperature thermal state is necessary because the energy gap between the singlet ground state and the threefold degenerate triplet first excited level practically closes around R_Li-H = 4Å, rendering the ground state alone an insufficient physical representation of the molecule. At lower separations (R_Li-H < 4 Å), the thermal state is effectively dominated by the ground state, as the gap remains nonzero. An extended analysis of entanglement in both the ground and excited states is provided in Supplementary Note 3.

Around equilibrium, there is still a considerable amount of entanglement between a Li-MEAO and a H-MEAO in the thermal state \(\widehat{\rho }(\beta )\), which closely resembles the ground state at this geometry. A fully ionic state would correspond to a product state with zero entanglement relative to the atomic partition. However, since the state retains some covalent character, the entanglement remains finite (see Supplementary Note 3 for the relationship between E_Li-H and the ionic character of the bond). Remarkably, the entanglement E_Li-H increases as the molecule dissociates, peaking around R_Li-H = 3 Å, coinciding with the avoided crossing. Beyond this peak, the entanglement decays to zero as dissociation progresses. This behavior stands in stark contrast to the monotonic decrease of bond order observed in Mulliken-like analyses⁴⁵, demonstrating how the MEAO formalism effectively captures subtle and complex bonding phenomena, such as the harpoon mechanism. These results underscore the strength of MEAOs as a comprehensive framework for analyzing challenging and nontrivial bonding scenarios.

Beyond Lewis structures: genuine multipartite entanglement

Genuine multipartite entanglement: While the majority of molecules conform to the Lewis paradigm of two-electron, two-center bonds, many molecules cannot be fully described by this model. These molecules often exhibit bonding structures involving more than two atomic centers. For two-center bonds, the connection between bonding and entanglement is established through the bipartite maximally entangled state \(\left|{\Psi }_{{{{\rm{bond}}}}}\right\rangle\) in Eq. (4). The natural question arises: what is the equivalent of \(\left|{\Psi }_{{{{\rm{bond}}}}}\right\rangle\) in the case of multicenter bonds?

To motivate our approach to the multicenter bonding problem, we first highlight that multipartite entangled states can belong to different classes with distinct internal structures. For example, for three qubits, the following two states⁴⁶:

belong to separate tripartite entanglement classes and are both maximally entangled within their respective classes. Notably, while every pair of qubits in \(\left|{{{\rm{W}}}}\right\rangle\) is entangled, no qubit pairs are entangled in \(\left|{{{\rm{GHZ}}}}\right\rangle\). Instead, the entanglement in \(\left|{{{\rm{GHZ}}}}\right\rangle\) exists collectively among all three qubits.

A K-partite pure state \(\left|\Psi \right\rangle\) is said to exhibit genuine multipartite entanglement (GME) if it cannot be expressed as a product state under any bipartition⁴⁷. Otherwise, the state is called biseparable. To determine whether \(\left|\Psi \right\rangle\) contains GME, one can check whether the entropy of any subsystem is nonzero. Based on this principle, a measure for pure state GME is defined as^48,49:

where A runs over all possible subsystems, which, in our case, correspond to collections of orbitals. When K = 2, GME reduces to the standard bipartite entanglement measure Eq. (5) for pure states. The GME measure Eq. (13) is normalized by its maximal value \(\log (d)\), where d is the dimension of the Fock space of the smallest subsystem. This measure can also be extended to mixed states via the convex roof construction^50,51,52. For example, the states \(\left|{{{\rm{GHZ}}}}\right\rangle\) and \(\left|{{{\rm{W}}}}\right\rangle\) have distinct GME values of \(\log (2)\) and \(\log (3)-\frac{2}{3}\log (2)\), respectively, demonstrating the capacity of GME to quantify entanglement genuinely shared among multiple orbitals, regardless of the specific multicenter bonding type or internal entanglement structure.

To illustrate the utility of GME as a multicenter bonding index, we analyze prototypical three-center, two- and four-electron bonds in the ethyl cation (\({{{{\rm{C}}}}}_{2}{{{{\rm{H}}}}}_{5}^{+}\)) and the allyl anion (\({{{{\rm{C}}}}}_{3}{{{{\rm{H}}}}}_{5}^{-}\)), respectively. In both molecules, a three-orbital cluster is detected from the MEAO correlation graph based on a B3LYP calculation, and the cluster is treated with a CAS. The remaining two-orbital, two-electron bonds are treated with mean-field accuracy, which ensures the correct electron count in the CAS. The MEAO framework thereby enables the automatic detection of electron-deficient and hypervalent bonds. As shown in Table 2, the constructed CAS directly identifies the three-center, two-electron bond in \({{{{\rm{C}}}}}_{2}{{{{\rm{H}}}}}_{5}^{+}\) and the three-center, four-electron bond in \({{{{\rm{C}}}}}_{3}{{{{\rm{H}}}}}_{5}^{-}\). The multicenter bonding character of these molecules are confirmed by the high GME values (0.891 for \({{{{\rm{C}}}}}_{2}{{{{\rm{H}}}}}_{5}^{+}\) and 0.787 for \({{{{\rm{C}}}}}_{3}{{{{\rm{H}}}}}_{5}^{-}\), both normalized by \(\log (4)\), as are all subsequent mentions of explicit values of GME). To further validate the universal applicability of GME as a multicenter bond index beyond carbon rings, we also analyzed lithium clusters. While lithium dimers are described by standard two-center two-electron bonds, lithium trimer anions (\({{{{\rm{Li}}}}}_{3}^{+}\)) and tetramers Li₄ display clear multicenter bonding patterns, with GME value 0.910 and 0.904, respectively.

GME as an aromaticity index: Aromaticity is a subtle chemical concept⁵³, manifested in many chemical and physical properties^54,55, but its underlying origin lies in the highly delocalized nature of electrons within a ring⁵⁶. This delocalization, a hallmark of aromatic systems, is naturally captured by the GME.

The most prominent example of an aromatic molecule is benzene, where the six out-of-plane p-orbitals collectively form a highly delocalized six-center bond, reflected by a six-π-orbital cluster with a GME value of 0.970 with the experimental geometry⁵⁷ and the cc-pVDZ basis set. The mutual information between the MEAOs of benzene is presented in Fig. 4, showing a distinct six-orbital cluster formed by six π-MEAOs. A complete active space configuration interaction calculation involving six electrons in six orbitals was performed for this cluster, while all other orbitals were treated at the Hartree-Fock level. In contrast to benzene, the highest six-orbital GME values for three six-member rings — cyclohexane (no π bond), cyclohexene (one π bond), and cyclohexa-1,3-diene (two π bonds) — are 0.015, 0.035, and 0.039, respectively. These results confirm that n-center electron delocalization is a necessary condition for an n-partite GME value to approach one. To benchmark GME as a quantitative aromaticity index, we selected challenging systems from a collection of aromaticity tests⁵⁸.

The alternative text for this image may have been generated using AI.

Full size image

Values of the normalized single orbital entropy \(S({\widehat{\rho }}_{i})/\log (4)\) and the normalized orbital-orbital correlation \({I}_{ij}/\log (16)\) are represented in log-scale by the color of the nodes and edges of the graph, respectively. The calculation is performed with the experimental equilibrium geometry and with the cc-pVDZ basis.

First, we analyze nitrogen-substituted benzene rings. Replacing one or more carbon atoms in the benzene ring with nitrogen disrupts the uniform electron delocalization, reducing aromaticity⁵⁹. This reduction is reflected in the GME values listed in Table 2, where the aromaticity of substituted species is slightly lower than that of benzene. The exact ordering of aromaticity for these six-membered rings differs depending on the aromaticity index used, such as the multicenter index (MCI)⁶⁰, harmonic oscillator measure of aromaticity (HOMA)^61,62, aromatic fluctuation index (FLU)⁶³, and nucleus-independent chemical shift (NICS)⁶⁴ (see Supplementary Note 4 for definitions). Notably, some well-known indices, such as HOMA and the para-delocalization index, failed this test⁵⁸. The GME values, however, correctly capture the trend that benzene has the highest aromaticity compared to nitrogen-substituted rings.

Second, we examine five-member rings, which can host five-orbital six-electron singlet bonds. A canonical example is \({{{{\rm{C}}}}}_{5}{{{{\rm{H}}}}}_{5}^{-}\), whose aromaticity is lower than benzene due to reduced point-group symmetry and the lack of particle-hole symmetry in the π-subspace. As expected, the GME value for \({{{{\rm{C}}}}}_{5}{{{{\rm{H}}}}}_{5}^{-}\) is lower than that for benzene. Substituting one carbon atom with nitrogen reduces the GME to 0.725, while substituting with oxygen reduces it further to 0.572. For comparison, the non-aromatic C₅H₆ yields a low GME value, and the detected cluster contains only four orbitals. These results confirm that GME accurately captures different levels of aromaticity in five-member rings.

Third, we investigate the effect of geometric deformation on the GME of benzene. The benzene ring is highly stable due to the resonance energy associated with its aromatic structure. Consequently, deviations from its equilibrium geometry are expected to raise the ground state energy and reduce aromaticity⁶⁵. We analyze six types of deformations: (1) bond length alternation (BLA), where C-C bond lengths alternate with a difference of ΔR; (2) bond length elongation (BLE) by ΔR; (3) clamping; (4) boat; (5) chair; and (6) pyramid (see Fig. 5a for graphical depictions). The GME values under these deformations are presented in Fig. 5b-e.

The alternative text for this image may have been generated using AI.

Full size image

a Six types of deformations of the benzene ring, including bond length alternation (BLA), bond length elongation (BLE), chair, boat, clamping, and pyramid. b GME as a function of the displacement ΔR in the BLA test. c GME as a function of the displacement ΔR in the BLE test. d GME as a function the deformation angle α in the chair (black circle), boat (blue square), and clamping (red triangle) test. e GME as a function of the deformation angle α in the pyramid test. f Energy (subtracted by the transition state energy, in atomic units, a.u.) of the ground state of the reaction systems along the intrinsic reaction coordinates (IRC) of a Diels-Alder cycloaddition. g Aromatic fluctuation index (FLU, grey circle), multicenter index (MCI, grey square), harmonic oscillator measure of aromaticity (HOMA, grey triangle), nucleus-independent chemical shift (NICS(1), grey diamond), and GME (red circle filled) of the ground state of the reaction systems along the IRC. Aromatic indices are linearly transformed to the interval [0, 1] for visibility. GME values are normalized by \(\log (4)\).

The GME is most sensitive to BLA, which explicitly breaks the six-fold symmetry of the C-C bond lengths. For BLE, we find that the GME decreases monotonically as the C-C distance increases, due to reduced overlap between the π-orbitals. However, the reduction in GME under BLE is less pronounced compared to BLA. Other deformations also reduce GME, as expected, with the exception of the pyramid deformation, where the GME slightly increases. This anomaly likely arises because the pyramid deformation preserves the six-fold symmetry of the ring. Overall, GME agrees well with the expectation that deviations from the equilibrium geometry reduce aromaticity, reinforcing its validity as a quantitative aromaticity index.

Aromaticity in transition states: During chemical reactions, electrons can become highly delocalized, sometimes forming transient aromatic ring structures. A concrete example is the Diels-Alder cycloaddition^66,67, where a conjugated butadiene reacts with ethene to form cyclohexene. Although neither the reactants nor the product is aromatic, the six π-electrons involved are temporarily shared across the entire ring as the π-bonds in butadiene and ethene break to form two new σ-bonds and a new π-bond, thereby promoting aromaticity^68,69 (see Fig. 5f for a graphical description of the reaction).

In Fig. 5f,g, we present the ground state energy of the system, calculated using the hybrid functional B3LYP⁷⁰ with the 6-311++G(d,p) basis set along the reaction pathway parameterized by the intrinsic reaction coordinate (IRC). We also include the GME values of the ground states based on corresponding active space calculations, alongside several aromaticity indices (shifted and renormalized to the interval [0, 1]; details in Methods). The multicenter bonding clusters and active spaces are determined using the graph partitioning technique described in Methods. Negative and positive IRC values indicate directions toward the reactants and products, respectively.

The Diels-Alder reaction presents a stringent test for aromaticity indices, as some widely used measures, such as HOMA^61,62 and FLU⁶³, fail to detect aromaticity in the transition state⁵⁸. Remarkably, GME exhibits a distinctive peak corresponding to a six-orbital cluster, with the location of the peak aligning well with other successful indices⁵⁸. Data points where GME = 0 correspond to regions where no prominent six-member clusters are detected. The successful detection of aromaticity in the transition state underscores the robustness of GME as an aromaticity index, even for systems far from equilibrium.

In summary, we have introduced GME as an index for multicenter bonding, particularly aromaticity, motivated by its high value in benzene. Through a series of tests, we confirmed that GME accurately assesses aromaticity across various aromatic molecules, including substituted and deformed benzene rings, while recovering well-known chemical trends both at equilibrium and in transition states. What sets GME apart is its conceptual foundation, rooted in the same formalism that characterizes two-center bonds⁷¹. Practically, the MEAO framework enables fully automatic identification of multicenter bonding clusters and their intensities, requiring no manual intervention. This makes GME a powerful tool for exploring complex bonding scenarios with minimal effort.

Bonding in transition metal complexes

Our final demonstration of the MEAO framework involves the chromium hexacarbonyl, Cr(CO)₆, a prototypical transition-metal complex. As an isolated molecule, CO exhibits a characteristic triple bond, comprising one σ and two π components. Upon coordination of the six CO ligands to the Cr center, the C-O bond order decreases to a value between a double and a triple bond, in accordance with the well-known Dewar-Chatt-Duncanson mechanism of σ-donation and π-back-donation^72,73. This reduction in covalency is also quantitatively shown by the diminished electron sharing between the C and O atoms in a subsequent study⁷⁴. In Table 3, we present the correlation and entanglement between the bonding MEAOs in both the free carbon monoxide and the chromium hexacarbonyl, correctly confirming the reduction in the covalent bond order between C and O. The MEAOs are optimized based on an IAO Hilbert space partition constructed from the entire atomic basis set (6-31G(d,p) basis set for carbon and oxygen, and Roos augmented double-zeta atomic natural orbital basis for chromium). Correlation and entanglement between the MEAOs are evaluated from the one-particle reduced density matrix (1RDM) obtained in a B3LYP density-functional-theory calculation.

Turning the fuzzy concept of chemical bonding into quantitative descriptors is essential for advancing our understanding of molecular properties and reaction mechanisms. In this work, we introduced an automated framework for versatile bonding analysis, combining the intuitive principles of valence bond theory with the rigorous tools of quantum information. At the core of this framework are the maximally entangled atomic orbitals (MEAOs), whose entanglement patterns naturally recover both Lewis (two-center) and beyond-Lewis (multicenter) bonding structures across a wide range of molecules. Moreover, our approach successfully captures challenging bonding scenarios that elude some widely adopted bonding descriptors, including the harpoon mechanism in LiH dissociation and aromaticity in the transition state of a Diels-Alder reaction.

Our framework offers a transformative perspective on the nature of chemical bonding. Established concepts such as electron sharing and spin pairing are seamlessly integrated into the entanglement between (or among) spatially localized orbitals. The chemical significance of MEAOs, which are automatically hybridized and intuitively meaningful, underscores a deep connection between the quantum mechanical essence of bonding and the entanglement between atomic subspaces. Furthermore, by leveraging insights from entanglement theory, our framework provides a unified characterization of two-center and multicenter bonds using the same tools. Notably, the LiH example demonstrates that the less commonly used Hilbert space atomic partition can be just as effective for bonding analysis as real-space partitioning, provided the correct descriptors, orbital entanglement in this case, are employed. This finding suggests a paradigm shift, where Hilbert space approaches could emerge as strong alternatives to real-space methods, offering advantages such as scalability and reduced sensitivity to integration errors.

Our work suggests several promising directions for future research. First, applying the MEAO formalism to monitor bond-breaking and bond-formation processes in molecular dynamics would be of practical interest. Second, while this study focused primarily on aromaticity as a multicenter bonding example, the efficacy of genuine multipartite entanglement in detecting other multicenter bonds, such as those in boranes⁷⁵, agostic bonds⁷⁶, and hydrogen bonds⁷⁷, warrants further exploration. Lastly, extending the treatment of multicenter bonds beyond a single active space description represents a compelling challenge. This could be addressed from a quantum information perspective, through the development of advanced multipartite entanglement measures for mixed states, or from a quantum chemistry perspective, by designing wave function methods capable of handling multiple active spaces simultaneously, while representing each multicenter bond by a pure state.

In conclusion, our framework bridges the gap between intuitive chemical bonding theories and quantitative quantum descriptors, offering a unified, scalable, and insightful approach for analyzing a wide array of bonding scenarios. By integrating concepts from quantum chemistry and quantum information theory, this work sets the foundation for a deeper understanding of chemical bonding and its role in molecular behavior.

Here we present a detailed account of the construction and usage of MEAOs for bonding analysis.

Step 1: Hilbert space partitioning. To introduce the MEAOs, we begin with an atomic partition \(\Pi={\{{{{{\mathcal{H}}}}}^{(m)}\}}_{m=1}^{M}\) of the one-particle Hilbert space \({{{\mathcal{H}}}}{=\bigoplus }_{m=1}^{M}{{{{\mathcal{H}}}}}^{(m)}\), where M is the number of atoms in the molecule. Typically, any procedure for such Hilbert space partitioning identifies both the subspaces \({{{{\mathcal{H}}}}}^{(m)}\) corresponding to individual atoms m and the localized orbitals \(\{{\phi }_{i}^{(m)}\}\) that form an orthonormal basis for \({{{{\mathcal{H}}}}}^{(m)}\). Common examples of these localized orbitals include IAOs²⁵ and Meta-Löwdin localized orbitals⁷⁸. In principle, each atomic subspace \({{{{\mathcal{H}}}}}^{(m)}\) can be spanned by infinitely many possible orbital bases \({{{{\mathcal{B}}}}}^{(m)}=\{{\phi }_{i}^{(m)}\}\), with any two such bases being related by an orbital rotation that preserves the subspace \({{{{\mathcal{H}}}}}^{(m)}\). When we perform orbital optimizations in later steps, we always consider orbital rotations that preserves the Hilbert space partition Π only.

Step 2: Constructing MEAOs. Guided by our analysis of the idealized covalent bond in Section “Entanglement in a single covalent bond”, the central idea of our approach is to identify the distinctive orbital basis \({{{{\mathcal{B}}}}}_{{{{\rm{MEAO}}}}}\) that maximizes the sum of inter-center orbital entanglements, \({\sum }_{i < j}E({\widehat{\rho }}_{ij})\), where the summation include orbital pairs (i, j) belonging to different atoms only. While direct computation of the entanglement \(E({\widehat{\rho }}_{ij})\) defined in Eq. (7) for all orbital pairs in larger systems can become computationally expensive, this approach offers valuable insight into the nature of chemical bonding by prioritizing the most significant inter-center correlations. Specifically, in the single covalent bond state, Eq. (4), both the particle numbers and the spin of the orthogonal atomic orbitals are perfectly correlated, as reflected in the maximal values of the coherent terms \(\langle 0,\uparrow \downarrow | {\widehat{\rho }}_{ij}| \uparrow \downarrow,0\rangle\) and \(\langle \uparrow,\downarrow | {\widehat{\rho }}_{ij}| \downarrow,\uparrow \rangle\), where \({\widehat{\rho }}_{ij}=\left|{\Psi }_{{{{\rm{bond}}}}}\right\rangle \left\langle {\Psi }_{{{{\rm{bond}}}}}\right|\) (see Supplementary Note 1 for a more detailed analysis). Inspired by this observation, we construct a proxy objective function for orbital entanglement that preserves the key features of the MEAOs in the perfect bonding state \(\left|{\Psi }_{{{{\rm{bond}}}}}\right\rangle\):

whose maximization for a given molecular wave function \(\left|\Psi \right\rangle\) determines the sought-after orbital basis \({{{{\mathcal{B}}}}}_{{{{\rm{MEAO}}}}}\) (see Fig. 1b). Here,

are elements of the 2RDM \( \Gamma \), and specifically \( {\Gamma }_{j \uparrow,j \downarrow }^{i \uparrow,i \downarrow }=\langle 0,\uparrow \downarrow | {\widehat{\rho }}_{ij}| \uparrow \downarrow,0 \rangle \) and \({\Gamma }_{j\uparrow,i\downarrow }^{i\uparrow,j\downarrow }=\langle \uparrow,\downarrow | {\widehat{\rho }}_{ij}| \downarrow,\uparrow \rangle\). The role of the proxy function F_MEAO is not to replicate the exact orbital entanglement for all two-orbital reduced density matrices, but rather to ensure that its maximum identifies an orbital basis closely aligned with the one that maximizes the true entanglement, a behavior supported by the analytical evidence discussed above and further validated by the numerical results. In most cases where the bonding pattern is simple (e.g., typical Lewis structures), the wave function \(\left|\Psi \right\rangle\) we use to compute the 2RDM (and correspondingly the objective function F_MEAO) can be a mean field wave function from a Hartree-Fock or a density functional theory calculation. In more challenging cases such as the dissociation of LiH, a correlated ground state (in our case a matrix product state, MPS) should be used.

Step 3: Grouping MEAOs into bonds. Bonds are identified by grouping MEAOs into strongly entangled pairs or clusters, which would require evaluating the entanglement between every inter-center pair of orbitals. This computationally expensive process can be efficiently approximated using the mutual information^{29,41,79,80,81,82}:

which is computationally inexpensive and serves as an upper bound to the entanglement, I_ij≥E_ij⁴¹. I_ij is also referred to as the total correlation, or simply correlation, between the orbitals i and j. To group MEAOs into bonds, we construct first a correlation graph: Two orbitals are considered connected in the correlation graph if their correlation I_ij exceeds a threshold η = 10% of the maximum value \({I}_{\max }=2\log (4)\). Similar to the maximal value of entanglement, \({I}_{\max }\) is also attained by the covalent bond state, Eq. (4). Introducing a threshold η > 0 reduces the number of orbitals to be analyzed by removing those that do not contribute meaningfully to bonding, while choosing η small enough ensures all relevant orbitals are retained. The resulting clusters reliably correspond to two-center or multicenter bonds, as illustrated in Fig. 1c.

Step 4: Apply quantum information tools to each identified bond. Once the orbital pairs and clusters responsible for chemical bonding are identified, quantum information tools are applied to characterize each bond in greater detail. These include measures of bipartite and genuinely multipartite entanglement, which provide deeper insight into the underlying bonding structure. To this end, we employ correlated wave functions to achieve an accurate description of the ground state. Two-center bonds are characterized through the entanglement between the two orbitals forming a pair after the graph partition, whereas multicenter bonds are quantified by the GME among all orbitals in a bonding cluster containing more than two centers. To compute the entanglement between two MEAOs in a bond, we directly evaluate the matrix elements of the corresponding two-orbital reduced density matrix using an MPS ansatz. These matrix elements are components of particle reduced density matrices up to the fourth order. This approach remains computationally efficient, as only a small subset of these matrix elements is required due to symmetry constraints.

Due to the difficulty of computing bipartite entanglement for mixed states, the GME can be efficiently computed only for pure states. Therefore, to evaluate the GME, we employ a complete active space (CAS) wave function in which only the orbitals participating in the multicenter bond are included in the active space, while the remaining orbitals are transformed into closed and virtual orbitals as described in the following. We first partition the one-particle Hilbert space into an active subspace (comprising the multicenter bonding orbitals) and a non-active subspace (the remainder). Let γ^non−active denote the restriction of the full spin-free one-particle reduced density matrix (1RDM), \({\gamma }_{ij}={\sum }_{\sigma=\uparrow,\downarrow }\langle \Psi | {f}_{i\sigma }^{{{\dagger}} }{f}_{j\sigma }| \Psi \rangle\), to the non-active subspace. Since the wave function on the non-active orbitals corresponds to almost fully occupied local bonding orbitals and almost empty local antibonding orbitals, the eigenvalues of γ^non−active are close to either 0 or 2. We assign the natural orbitals of γ^non−active with eigenvalues greater than 1 as closed orbitals, and those with eigenvalues smaller than 1 as virtual orbitals. By construction, the closed and virtual orbitals are orthogonal to those involved in the multicenter bond, thereby defining a consistent active space partition.

Numerical details: The construction of intrinsic atomic orbitals and meta-Löwdin orbitals was performed using the PySCF package⁸³ with default settings. To find the maximally entangled atomic orbitals (MEAOs), the objective function, Eq. (14), was maximized using a second-order Newton-Raphson method, where the analytic gradient and diagonal elements of the Hessian were used, as detailed in Supplementary Note 1.

Molecular structures of the molecules in Table 1 and CO in Table 3 were taken from experimental data in the NIST Computational Chemistry Comparison and Benchmark Database⁵⁷, except for He₂ whose bond length is taken from ref. ⁸⁴. The carbon rings in Table 2 and those in Fig. 5b-e are optimized using the B3LYP functional with the 6-311++G(d,p) basis⁵⁸. In particular, the geometry optimization is performed with the PySCF package for \({{{{\rm{C}}}}}_{2}{{{{\rm{H}}}}}_{5}^{+}\), \({{{{\rm{C}}}}}_{3}{{{{\rm{H}}}}}_{5}^{+}\), C₆H₈, C₆H₁₀, and C₆H₁₂, while the optimized geometries of the rest of the carbon rings are directly taken from the Supplemental Material of ref. ⁵⁸. The geometries of \({{{{\rm{Li}}}}}_{3}^{+}\) and Li₄ are optimized geometries on CCSD and CCSDT levels of theory, respectively⁸⁵. The structure of Cr(CO)₆ is constructed with the experimental bond lengths reported in ref. ⁸⁶. All active space and full configuration interaction calculations were performed with an MPS and a density matrix renormalization group (DMRG) solver provided by Block2^87,88.

The relative entropy of bipartite entanglement for mixed states was calculated using the semidefinite programming package CVX^89,90 and quantum information helper functions provided in refs. ^91,92. The algorithm for calculating the entanglement can be found in ref. ⁹³.

The Diels-Alder reaction was simulated using the Gaussian program⁹⁴, starting from the transition state and progressing towards the reagents and products along separate paths. For the thermal state of LiH, the thermal 2RDM is used to optimize the MEAOs, starting from a set of meta-Löwdin localized orbitals. The thermal 2RDM is computed as the thermal average of the state-specific 2RDMs obtained via DMRG calculations correlating all electrons and all orbitals in the basis.

The minimization in the definition of GME is approximated by minimizing over all single orbital entropies, two orbital entropies, and bipartite entanglement along the MPS. This is a sensible approximation since the largest genuinely entanled cluster we have studied in this work consists of only six centers. In most cases, the minimization is realized by one of the single orbital entropy values. The aromatic indices in Fig. 5g were taken from ref. ⁵⁸ and linearly renormalized to the interval [0,1], by first taking the absolute value, then subtracting by the lowest value, and finally divided by the largest value. Orbital isosurface plots were produced using the software Jmol⁹⁵.

1. Lewis, G. N. The atom and the molecule. J. Am. Chem. Soc. 38, 762 (1916).

   Article  ADS  CAS  Google Scholar 

2. Heitler, W. & London, F. Wechselwirkung neutraler Atome und homöopolare Bindung nach der Quantenmechanik. Z. Phys. 44, 455 (1927).

   Article  ADS  CAS  Google Scholar 

3. Gerratt, J., Cooper, D., Karadakov, P. A. & Raimondi, M. Modern valence bond theory. Chem. Soc. Rev. 26, 87 (1997).

   Article  CAS  Google Scholar 

4. Shaik, S. S., & Hiberty, P. C. A Chemist’s Guide to Valence Bond Theory (John Wiley & Sons, 2007) https://doi.org/10.1002/9780470192597.

5. Roothaan, C. C. J. New developments in molecular orbital theory. Rev. Mod. Phys. 23, 69 (1951).

   Article  ADS  CAS  Google Scholar 

6. Malmqvist, P. -Å & Roos, B. O. The CASSCF state interaction method. Chem. Phys. Lett. 155, 189 (1989).

   Article  ADS  CAS  Google Scholar 

7. Andersson, K., Malmqvist, P. A., Roos, B. O., Sadlej, A. J. & Wolinski, K. Second-order perturbation theory with a CASSCF reference function. J. Phys. Chem. 94, 5483 (1990).

   Article  ADS  CAS  Google Scholar 

8. Schollwöck, U. The density-matrix renormalization group. Rev. Mod. Phys. 77, 259 (2005).

   Article  ADS  MathSciNet  Google Scholar 

9. Chan, G. K.-L. & Sharma, S. The density matrix renormalization group in quantum chemistry. Annu. Rev. Phys. Chem. 62, 465 (2011).

   Article  ADS  CAS  PubMed  Google Scholar 

10. Bartlett, R. & Musiał, M. Coupled-cluster theory in quantum chemistry. Rev. Mod. Phys. 79, 291 (2007).

    Article  ADS  CAS  Google Scholar 

11. Mulliken, R. S. Electronic population analysis on LCAO–MO molecular wave functions. I. J. Chem. Phys. 23, 1833 (1955).

12. Bader, R. F. & Stephens, M. E. Spatial localization of the electronic pair and number distributions in molecules. J. Am. Chem. Soc. 97, 7391 (1975).

    Article  ADS  CAS  Google Scholar 

13. Bridgeman, A. J., Cavigliasso, G., Ireland, L. R. & Rothery, J. The Mayer bond order as a tool in inorganic chemistry. Dalton Trans. 2095 https://doi.org/10.1039/B102094N (2001).

14. Mayer, I. Bond order and valence indices: a personal account. J. Comput. Chem. 28, 204 (2007).

    Article  ADS  CAS  PubMed  Google Scholar 

15. Roos, B. O., Borin, A. C. & Gagliardi, L. Reaching the maximum multiplicity of the covalent chemical bond. Angew. Chem. Int. Ed. 46, 1469 (2007).

    Article  CAS  Google Scholar 

16. Bader, R. F. Atoms in molecules. Acc. Chem. Res. 18, 9 (1985).

    Article  CAS  Google Scholar 

17. Silvi, B. & Savin, A. Classification of chemical bonds based on topological analysis of electron localization function. Nature 371, 683 (1994).

    Article  ADS  CAS  Google Scholar 

18. Becke, A. D. & Edgecombe, K. E. A simple mesure of electron localization in atomic and molecular systems. J. Comp. Phys. 92, 5397 (1990).

    CAS  Google Scholar 

19. Pendás, Á. M., & Contreras-García, J. Topological Approaches to the Chemical Bond (Springer Nature, 2023) https://doi.org/10.1007/978-3-031-13666-5.

20. Fradera, X., Austen, M. A. & Bader, R. F. The Lewis model and beyond. J. Phys. Chem. A 103, 304 (1999).

    Article  CAS  Google Scholar 

21. Matito, E. An electronic aromaticity index for large rings. Phys. Chem. Chem. Phys. 18, 11839 (2016).

    Article  CAS  PubMed  Google Scholar 

22. Van Hende, D., Van der Stichelen, R., Bultinck, P. & Acke, G. Quantum information patterns between atoms in a molecule. Chem. Eur. J. 30, e202400819 (2024).

    Article  PubMed  Google Scholar 

23. Mayer, I. Towards a “Chemical” Hamiltonian. Int. J. Quantum Chem. 23, 341 (1983).

    Article  CAS  Google Scholar 

24. Reed, A. E., Weinstock, R. B. & Weinhold, F. Natural population analysis. J. Chem. Phys. 83, 735 (1985).

    Article  ADS  CAS  Google Scholar 

25. Knizia, G. Intrinsic atomic orbitals: an unbiased bridge between quantum theory and chemical concepts. J. Chem. Theory Comput. 9, 4834 (2013).

    Article  CAS  PubMed  Google Scholar 

26. Mayer, I., Bakó, I. & Stirling, A. Are there atomic orbitals in a molecule? J. Phys. Chem. A 115, 12733 (2011).

    Article  CAS  PubMed  Google Scholar 

27. Nalewajski, R. F. & Parr, R. G. Information theory, atoms in molecules, and molecular similarity. Proc. Natl. Acad. Sci. USA. 97, 8879 (2000).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

28. Nalewajski, R. F. Information Theory of Molecular Systems (Elsevier, 2006) https://doi.org/10.1016/b978-0-444-51966-5.x5063-4.

29. Rissler, J., Noack, R. M. & White, S. R. Measuring orbital interaction using quantum information theory. Chem. Phys. 323, 519 (2006).

    Article  CAS  Google Scholar 

30. Boguslawski, K., Tecmer, P., Barcza, G., Legeza, O. & Reiher, M. Orbital entanglement in bond-formation processes. J. Chem. Theory Comput. 9, 2959 (2013).

    Article  CAS  PubMed  Google Scholar 

31. Krumnow, C., Veis, L., Legeza, Ö & Eisert, J. Fermionic orbital optimization in tensor network states. Phys. Rev. Lett. 117, 210402 (2016).

    Article  ADS  CAS  PubMed  Google Scholar 

32. Stein, C. J. & Reiher, M. Automated selection of active orbital spaces. J. Chem. Theory Comput. 12, 1760 (2016).

    Article  CAS  PubMed  Google Scholar 

33. Stein, C. J., Pantazis, D. A. & Krewald, V. Orbital entanglement analysis of exchange-coupled systems. J. Phys. Chem. Lett. 10, 6762 (2019).

    Article  CAS  PubMed  Google Scholar 

34. Ding, L., Knecht, S., Zimborás, Z. & Schilling, C. Quantum correlations in molecules: from quantum resourcing to chemical bonding. Quantum Sci. Technol. 8, 015015 (2022).

    Article  ADS  Google Scholar 

35. Bensberg, M. & Reiher, M. Corresponding active orbital spaces along chemical reaction paths. J. Phys. Chem. Lett. 14, 2112 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

36. Ding, L., Knecht, S. & Schilling, C. Quantum information-assisted complete active space optimization (QICAS). J. Phys. Chem. Lett. 14, 11022 (2023).

    Article  CAS  PubMed  Google Scholar 

37. Liao, K., Ding, L. & Schilling, C. Quantum information orbitals (QIO): unveiling intrinsic many-body complexity by compressing single-body triviality. J. Phys. Chem. Lett. 15, 6782 (2024).

38. Cotton, F. A. Chemical Applications of Group Theory (John Wiley & Sons, 1991).

39. Vedral, V., Plenio, M. B., Rippin, M. A. & Knight, P. L. Quantifying entanglement. Phys. Rev. Lett. 78, 2275 (1997).

    Article  ADS  MathSciNet  CAS  Google Scholar 

40. Henderson, L. & Vedral, V. Information, relative entropy of entanglement, and irreversibility. Phys. Rev. Lett. 84, 2263 (2000).

    Article  ADS  CAS  PubMed  Google Scholar 

41. Modi, K., Paterek, T., Son, W., Vedral, V. & Williamson, M. Unified view of quantum and classical correlations. Phys. Rev. Lett. 104, 080501 (2010).

    Article  ADS  MathSciNet  PubMed  Google Scholar 

42. Matito, E., Solà, M., Salvador, P. & Duran, M. Electron sharing indexes at the correlated level. application to aromaticity calculations. Faraday Discuss. 135, 325 (2007).

    Article  ADS  CAS  PubMed  Google Scholar 

43. Polanyi, J. C. & Zewail, A. H. Direct observation of the transition state. Acc. Chem. Res. 28, 119 (1995).

    Article  CAS  Google Scholar 

44. Rodríguez-Mayorga, M., Ramos-Cordoba, E., Salvador, P., Solà, M. & Matito, E. Bonding description of the harpoon mechanism. Mol. Phys. 114, 1345 (2016).

    Article  ADS  Google Scholar 

45. Ponec, R. & Cooper, D. L. Anatomy of bond formation. bond length dependence of the extent of electron sharing in chemical bonds. J. Mol. Struct. THEOCHEM 727, 133 (2005).

    Article  CAS  Google Scholar 

46. Dür, W., Vidal, G. & Cirac, J. I. Three qubits can be entangled in two inequivalent ways. Phys. Rev. A 62, 062314 (2000).

    Article  ADS  MathSciNet  Google Scholar 

47. Horodecki, M., Horodecki, P. & Horodecki, R. Separability of n-particle mixed states: necessary and sufficient conditions in terms of linear maps. Phys. Lett. A 283, 1 (2001).

    Article  ADS  MathSciNet  CAS  Google Scholar 

48. Pope, D. T. & Milburn, G. J. Multipartite entanglement and quantum state exchange. Phys. Rev. A 67, 052107 (2003).

    Article  ADS  Google Scholar 

49. Ma, Z.-H. et al. Measure of genuine multipartite entanglement with computable lower bounds. Phys. Rev. A 83, 062325 (2011).

    Article  ADS  Google Scholar 

50. Wootters, W. K. Entanglement of formation of an arbitrary state of two qubits. Phys. Rev. Lett. 80, 2245 (1998).

    Article  ADS  CAS  Google Scholar 

51. Nielsen, M. A. Continuity bounds for entanglement. Phys. Rev. A 61, 064301 (2000).

    Article  ADS  MathSciNet  Google Scholar 

52. Wei, T.-C. & Goldbart, P. M. Geometric measure of entanglement and applications to bipartite and multipartite quantum states. Phys. Rev. A 68, 042307 (2003).

    Article  ADS  Google Scholar 

53. Solà, M. Why aromaticity is a suspicious concept? Why? Front. Chem. 5, 22 (2017).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

54. Schleyer, P. V. R. & Jiao, H. What is aromaticity. Pure Appl. Chem. 68, 209 (1996).

    Article  ADS  CAS  Google Scholar 

55. Schleyer, P. V. R. Introduction: aromaticity. Chem. Rev. 101, 1115 (2001).

    Article  CAS  PubMed  Google Scholar 

56. Feixas, F., Matito, E., Poater, J. & Solà, M. Quantifying aromaticity with electron delocalisation measures. Chem. Soc. Rev. 44, 6389 (2015).

    Article  Google Scholar 

57. NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, http://cccbdb.nist.gov/,release 22, May 2022, Editor: Russell D. Johnson III.

58. Feixas, F., Matito, E., Poater, J. & Solà, M. On the performance of some aromaticity indices: a critical assessment using a test set. J. Comp. Chem. 29, 1543 (2008).

    Article  ADS  CAS  Google Scholar 

59. Hoffmann, R. Extended Hückel Theory. II. σ Orbitals in the Azines. J. Chem. Phys. 40, 2745 (1964).

    Article  ADS  CAS  Google Scholar 

60. Bultinck, P. et al. Electron delocalization and aromaticity in linear polyacenes: atoms in molecules multicenter delocalization index. J. Phys. Chem. A 110, 7642 (2006).

    Article  CAS  PubMed  Google Scholar 

61. Kruszewski, J. & Krygowski, T. Definition of aromaticity basing on the harmonic oscillator model. Tetrahedron Lett. 13, 3839 (1972).

    Article  Google Scholar 

62. Krygowski, T. M. Crystallographic studies of inter-and intramolecular interactions reflected in aromatic character of π-electron systems. J. Chem. Inf. Comput. Sci. 33, 70 (1993).

63. Matito, E., Duran, M. & Sola, M. The aromatic fluctuation index (FLU): A new aromaticity index based on electron delocalization. J. Chem. Phys. 122, 014109 (2005).

64. Schleyer, P.vR., Maerker, C., Dransfeld, A., Jiao, H. & van Eikema Hommes, N. J. Nucleus-independent chemical shifts: a simple and efficient aromaticity probe. J. Am. Chem. Soc. 118, 6317 (1996).

    Article  ADS  CAS  PubMed  Google Scholar 

65. Feixas, F., Matito, E., Poater, J. & Sola, M. Aromaticity of distorted benzene rings: exploring the validity of different indicators of aromaticity. J. Phys. Chem. A 111, 4513 (2007).

    Article  CAS  PubMed  Google Scholar 

66. Diels, O. & Alder, K. Synthesen in der hydroaromatischen Reihe. Justus Liebigs Ann. Chem. 460, 98 (1928).

    Article  CAS  Google Scholar 

67. Nicolaou, K. C., Snyder, S. A., Montagnon, T. & Vassilikogiannakis, G. The Diels–Alder reaction in total synthesis. Angew. Chem. Int. Ed. 41, 1668 (2002).

    Article  CAS  Google Scholar 

68. Manoharan, M., De Proft, F. & Geerlings, P. A computational study of aromaticity-controlled Diels–Alder reactions. J. Chem. Soc., Perkin Trans. 2, 1767–1773 (2000).

69. Matito, E., Poater, J., Duran, M. & Solà, M. An analysis of the changes in aromaticity and planarity along the reaction path of the simplest Diels–Alder reaction. Exploring the validity of different indicators of aromaticity. J. Mol. Struct. THEOCHEM 727, 165 (2005).

    Article  CAS  Google Scholar 

70. Stephens, P. J., Devlin, F. J., Chabalowski, C. F. & Frisch, M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. J. Phys. Chem. 98, 11623 (1994).

    Article  ADS  CAS  Google Scholar 

71. Ponec, R. & Uhlik, F. Multicentre bond indices from the generalized population analysis of higher order densities. Croat. Chem. Acta 69, 941 (1996).

    CAS  Google Scholar 

72. Dewar, J. A review of the π-complex theory. Bull. Soc. chim. Fr. 18, C71 (1951).

    Google Scholar 

73. Chatt, J. & Duncanson, L. A. 586. Olefin co-ordination compounds. Part III. Infra-red spectra and structure: attempted preparation of acetylene complexes. J. Chem. Soc. 2939 https://pubs.rsc.org/en/content/articlelanding/1953/jr/jr9530002939 (1953).

74. Matito, E. & Solà, M. The role of electronic delocalization in transition metal complexes from the electron localization function and the quantum theory of atoms in molecules viewpoints. Coord. Chem. Rev. 253, 647 (2009).

    Article  CAS  Google Scholar 

75. Melichar, P., Hnyk, D. & Fanfrlík, J. A systematic examination of classical and multi-center bonding in heteroborane clusters. Phys. Chem. Chem. Phys. 20, 4666 (2018).

    Article  CAS  PubMed  Google Scholar 

76. Brookhart, M. & Green, M. L. Carbon hydrogen-transition metal bonds. J. Organomet. Chem. 250, 395 (1983).

    Article  CAS  Google Scholar 

77. Giambiagi, M., de Giambiagi, M. S. & Mundim, K. C. Definition of a multicenter bond index. Struct. Chem. 1, 423 (1990).

    Article  CAS  Google Scholar 

78. Sun, Q. & Chan, G. K.-L. Exact and optimal quantum mechanics/molecular mechanics boundaries. J. Chem. Theory Comput. 10, 3784–3790 (2014).

79. Adami, C. & Cerf, N. J. von Neumann capacity of noisy quantum channels. Phys. Rev. A 56, 3470 (1997).

    Article  ADS  CAS  Google Scholar 

80. Groisman, B., Popescu, S. & Winter, A. Quantum, classical, and total amount of correlations in a quantum state. Phys. Rev. A 72, 032317 (2005).

    Article  ADS  MathSciNet  Google Scholar 

81. Boguslawski, K. & Tecmer, P. Orbital entanglement in quantum chemistry. Int. J. Quantum Chem. 115, 1289 (2015).

    Article  CAS  Google Scholar 

82. Ding, L. et al. Concept of orbital entanglement and correlation in quantum chemistry. J. Chem. Theory Comput. 17, 79 (2020).

    Article  PubMed  Google Scholar 

83. Sun, Q. et al. Pyscf: the python-based simulations of chemistry framework. Wiley Interdiscip. Rev. Comput. Mol. Sci. 8, e1340 (2018).

    Article  Google Scholar 

84. Grisenti, R. E. et al. Determination of the bond length and binding energy of the helium dimer by diffraction from a transmission grating. Phys. Rev. Lett. 85, 2284 (2000).

    Article  ADS  CAS  PubMed  Google Scholar 

85. Rousseau, R. & Marx, D. Ab initio calculations on small lithium clusters. Phys. Rev. A 56, 617 (1997).

    Article  ADS  CAS  Google Scholar 

86. Jost, A., Rees, B. & Yelon, W. Electronic structure of chromium hexacarbonyl at 78 KI Neutron diffraction study. Acta Crystallogr. B. 31, 2649 (1975).

    Article  ADS  Google Scholar 

87. Zhai, H. & Chan, G. K. Low communication high performance ab initio density matrix renormalization group algorithms. J. Chem. Phys. 154, 224116 (2021).

    Article  ADS  CAS  PubMed  Google Scholar 

88. Zhai, H. et al. Block2: A comprehensive open source framework to develop and apply state-of-the-art DMRG algorithms in electronic structure and beyond. J. Chem. Phys. 159, 234801 (2023).

89. Grant, M. & Boyd, S. CVX: Matlab software for disciplined convex programming, version 2.1, https://cvxr.com/cvx (2014).

90. Grant, M. & Boyd, S. Graph implementations for nonsmooth convex programs. in Recent Advances in Learning and Control, Lecture Notes in Control and Information Sciences, edited by V. Blondel, S. Boyd and H. Kimura (Springer-Verlag Limited, 2008) pp. 95–110,http://stanford.edu/boyd/graph_dcp.html.

91. Johnston, N. QETLAB: A MATLAB toolbox for quantum entanglement, version 1.0, https://qetlab.com (2016).

92. Fawzi, H., Saunderson, J. & Parrilo, P. A. Semidefinite approximations of the matrix logarithm. Foundations of Computational Mathematics (2018) package cvxquad at https://github.com/hfawzi/cvxquad.

93. Ding, L. Fermionic entanglement and correlation arXiv:2207.03848 https://arxiv.org/abs/2207.03848 (2022).

94. Frisch, M. J. et al. Gaussian 16 Revision B.01 (2016), Gaussian Inc. Wallingford CT.

95. Jmol: an open-source Java viewer for chemical structures in 3D.Stop http://www.jmol.org/.

96. Ding, L., Matito, E., & Schilling, C. Maximally entangled atomic orbitals, (2026) https://doi.org/10.5281/zenodo.19553976.

Download references

Authors and Affiliations

1. Faculty of Physics, Arnold Sommerfeld Centre for Theoretical Physics (ASC), München, Germany

   Lexin Ding & Christian Schilling

2. Munich Center for Quantum Science and Technology (MCQST), München, Germany

   Lexin Ding & Christian Schilling

3. Donostia International Physics Center (DIPC), Donostia, Spain

   Eduard Matito

4. Ikerbasque Foundation for Science, Bilbao, Spain

   Eduard Matito

5. Donostia Polimero eta Material Aurreratuak: Fisika, Kimika eta Teknologia, Kimika Fakultatea, Euskal Herriko Unibertsitatea (EHU), Donostia, Spain

   Eduard Matito

Contributions

The project was conceived by C.S. The theoretical framework was developed by L.D. with input from C.S., and the codes were implemented by L.D.; E.M. provided expertise in chemical bonding theory and aromaticity, designed most of the examples in the manuscript and contributed to the chemical interpretation of the results. The manuscript was drafted by L.D. and revised by all authors (L.D., E.M., and C.S.). Funding was acquired by C.S.

Corresponding author

Correspondence to Christian Schilling.

Competing interests

The authors declare no competing interests.

Peer review information

Nature Communications thanks the anonymous reviewers for their contribution to the peer review of this work. A peer review file is available

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions