The location of the chemical bond. Application of long covalent bond theory to the structure of silica

Oxygen is the most abundant terrestrial element and is found in a variety of materials, but still wanting is a universal theory for the stability and structural organization it confers. Herein, a computational molecular orbital analysis elucidates the structure, stability, and cooperative bonding of α-quartz silica (SiO2). Despite geminal oxygen-oxygen distances of 2.61–2.64 Å, silica model complexes exhibit anomalously large O-O bond orders (Mulliken, Wiberg, Mayer) that increase with increasing cluster size—as the silicon-oxygen bond orders decrease. The average O-O bond order in bulk silica computes to 0.47 while that for Si-O computes to 0.64. Thereby, for each silicate tetrahedron, the six O-O bonds employ 52% (5.61 electrons) of the valence electrons, while the four Si-O bonds employ 48% (5.12 electrons), rendering the O-O bond the most abundant bond in the Earth’s crust. The isodesmic deconstruction of silica clusters reveals cooperative O-O bonding with an O-O bond dissociation energy of 4.4 kcal/mol. These unorthodox, long covalent bonds are rationalized by an excess of O 2p–O 2p bonding versus anti-bonding interactions within the valence molecular orbitals of the SiO4 unit (48 vs. 24) and the Si6O6 ring (90 vs. 18). Within quartz silica, oxygen 2p orbitals contort and organize to avoid molecular orbital nodes, inducing the chirality of silica and resulting in Möbius aromatic Si6O6 rings, the most prevalent form of aromaticity on Earth. This long covalent bond theory (LCBT) relocates one-third of Earth’s valence electrons and indicates that non-canonical O-O bonds play a subtle, but crucial role in the structure and stability of Earth’s most abundant material.


Introduction
Silica (SiO 2 ) constitutes about 59% of the Earth's crust (Clarke and Washington, 1924) and its abundance is largely explained by its stability. By mass, silica is 53% oxygen, the most abundant element on Earth. A universal connection between oxygen and the material stability it bestows has not been described. For years, Pauling argued that a "molecule is stabilized by. . . resonance, its energy being less than the energy corresponding to any one of the structures among which it resonates." (Pauling, 1946) "It is this extra stability of the system. . . that is called the resonance energy." (Pauling, 1960) Although Pauling never quantified this value for silica (Torgunrud et al., 2020) his detailed resonance description of "low quartz", illustrated in Figure 1, distinctly implies that it enjoys a resonance energy (Pauling, 1980). The 50% ionic/50% covalent character of silica was ardently defended by Pauling and determined by the electronegativity difference of Δx = 1.7 between oxygen (x = 3.5) and silicon (x = 1.8) (Pauling, 1952). With the premise of 50% ionic character and only siliconoxygen single bonds, resonance hybrid averaging cannot achieve the charges of +1 on silicon and −0.5 on each oxygen, as limited by the electroneutrality postulate (Newton and Gibbs, 1980). However, Pauling's silicon-oxygen double bonds introduce negative charge on silicon and positive charge on oxygen, allowing for resonance averaging to Si +1 and O −0.5 when the ionic character is 50%. Since each silicon atom forms two single bonds and two double bonds, Pauling concluded that the Si-O bond order in silica is 3/2, with the double bonds accommodated by sp 3 d 2 orbitals of silicon (Pauling, 1980). While a topic of considerable debate (Janes and Oldfield, 1986), recent experimental and computational results indicate that O 2p-Si 3d π-bonding is minimal in silica (Silvi et al., 1997;Garvie et al., 2000;Gibbs et al., 2009) and that Si-O bond orders are less than unity (Newton and Gibbs, 1980;Noritake and Kawamura, 2015;Torgunrud et al., 2020), contradicting Pauling's resonance scheme and stated bond order. Without silicon-oxygen double bonds, the 50% covalent component defaults to the single silicon-oxygen bonds of the canonical resonance form. Resonance averaging that is 50% canonical and 50% ionic yields a net +2 charge on silicon and an obvious violation of the electroneutrality postulate.
Herein, an alternative resonance formulation for silica is presented that ameliorates these incongruities. Non-canonical, long covalent oxygen-oxygen bonds are proposed and rationalized-ostensibly necessitated-by a molecular orbital analysis of silica model complexes, thereby demystifying the structure, bonding, and stability of Earth's most abundant material.
2 Results and discussion

Resonance hybrids with long bonds
The ionic component of Pauling's 1980 resonance scheme for silica ( Figure 2A) employs oxidized silicon(IV) and reduced oxygen(-II); the covalent component employs reduced silicon(-II) and oxidized oxygen(I). An alternative covalent resonance formulation is proposed in Figure 2B, which assumes the same oxidation states, but a different arrangement of bonds. Instead of arguably discounted silicon-oxygen double bonds, long oxygen-oxygen bonds averaging 2.63 Å are enlisted. Each silicon-oxygen interaction comprises a 50% single bond and 50% no-bond (Roberts et al., 1950) resonance formulation; hence the predicted Si-O bond order is 1/2. Each geminal oxygen-oxygen interaction comprises a 33% single bond and 67% no-bond resonance formulation; hence the predicted O-O bond order is 1/3. Note that this resonance scheme involves an equal number of Si-O (4 = 4 x 2 x ½) and O-O (4 = 6 x 2 x ⅓) valence electrons.

Bond order analysis
This non-canonical O-O bond can be interrogated by a Mulliken bond order analysis (Mulliken, 1955a), which "characterizes the accumulation of the electrons in the region between the chemically bonded atoms, and is a very useful quantity often characterizing well the bond strength" (Mayer, 2007) and bond covalency (Mulliken, 1955b;Segall et al., 1996). The canonical silica structure employs 8 of its 16 valence electrons (= 4 + 6 + 6) to make four Si-O single bonds, leaving 4 lone pairs of non-bonding electrons on the oxygen atoms. Computed at the DFT/B3LYP/6-311++G** level of theory (Shao et al., 2015;Levien et al., 1980) 1 , silicic acid (Si(OH) 4 ) appears nearly canonical with four Si-O bond orders averaging 1.03 (range = 1.01-1.04) and six oxygen-oxygen bond orders averaging an inconspicuous 0.05 (range = 0.04-0.06). Accordingly, each silicic acid drafts 8.21 valence electrons for Si-O bonding and 0.55 electrons for O-O bonding, totaling 8.77 valence bonding electrons and leaving 7.23 as non-bonding.
However, for the silicate at the center of a large 29-silicon cluster (Si29, Si 29 O 76 H 36 ), where structural regularity is high and edge effects are low, the Si-O bond orders are substantially less than unity at 0.63, 0.63, 0.68, and 0.68, averaging 0.66 and requiring only 5.25 electrons. The O-O bond orders are conspicuously greater than zero at 0.38, 0.38, 0.43, 0.43, 0.63, and 0.63, averaging 0.48 and requiring 5.78 electrons-with these six O-O bonds occupying the edges of a silicate tetrahedron, as shown in Figure 3A. Thus, the central silicate of the Si29 cluster employs 11.03 valence bonding electrons, leaving 4.97 electrons as non-bonding. Moreover, this central silicate enlists 5.23 (= 5.78-0.55) additional electrons for oxygen-oxygen bonding versus silicic acid. The electrostatic charge on silicon is +0.94e and the charges on the oxygen atoms are −0.51, −0.51, −0.54, and −0.54e, averaging −0.52e. The computed bond orders and atomic charges mirror the simplistic resonance hybrid averaging of Figure 2B. Figure 3B plots the computed central silicate Si-O and O-O bond orders versus the number of silicon atoms in model quartz silica clusters ranging from Si1 (silicic acid) to Si35. As the model complex is enlarged from Si1 to Si5, the average O-O bond order increases about 0.05 for each additional silicate unit appended to the central silicate. Thereafter, the O-O bond order increases about 0.01 for each additional silicate-reaching 0.48 for Si29 and 0.47 for Si35 (presently the computational limit). The average Si-O bond order decreases with increasing cluster size-albeit with less regularity-reaching 0.66 for Si29 and 0.59 for Si35. Note also that the Si 3d atomic orbital contribution to the valence molecular orbitals is consistently below 1% (Si1, 0.8%; Si5, 0.7%; Si29, 0.4%), further bolstering the premise that Si-O π bonding is minimal in silica. For these same model clusters, Figure 3C plots the percentage of valence electrons belonging to the Si-O and O-O bonds of the central silicate unit. The O-O valence electron percentage climbs steeply to 36% (Si5) as the first shell of four silicates is added to the central silicate. The second shell accommodates twelve additional silicates, but their addition has a less dramatic effect. With 2-7 second-shell silicates (Si8-Si18), the O-O valence electron percentage spans 37%-43%. With 8-12 second-shell silicates (Si21-Si35), the percentage spans 51%-54%. These plots affirm a

FIGURE 2
An alternative to Pauling's resonance scheme for silica (A) proposes long oxygen-oxygen bonds, ranging from 2.61 to 2.64 Å (B). This alternative conforms to Pauling's 50% ionic/50% covalent apportionment, enlists the same average atomic charges and oxidation states, avoids unwarranted  (Kirfel and Eichhorn, 1990), the second most abundant material in the Earth's crust (15%) (Clarke and Washington, 1924

Molecular orbital analysis
Generally, atoms are not covalently bonded unless they experience an excess of bonding interactions versus anti-bonding interactions and thereby, "a concentration of charge between the two nuclei." (Feynman, 1939) Similarly, Hoffmann stated succinctly that "positive overlap implies stabilization or bonding." (Hoffmann, 1971) So, is there charge concentration and positive overlap between oxygen atoms of silica that account for the non-canonical O-O bond orders of Figure 3B?
The repeat unit of silica can be considered as SiO 2 or Si(½O) 4 and the four oxygen atoms of the latter can be modeled by Si(OSi) 4 to provide the molecular orbitals (MOs) of Figure 4A. Depicted are sixteen valence MOs built primarily from the four oxygen 2s orbitals and the twelve oxygen 2p orbitals (IsoValue = 0.01 [e/bohr 3 ] 0.5 ). The oxygen character is above 70% for eleven of these MOs, averages 71%, and exceeds the valence electron percentage contributed by oxygen, 55% 2 . The remaining six valence MOs (HOMO-5 through HOMO, not shown) are built primarily from silicon atomic orbitals and have a silicon character averaging 86%-well exceeding the valence electron percentage contributed by silicon, 45%. Clearly, there is disproportional and poor silicon/oxygen atomic orbital mixing and the oxygen-rich MOs rationalize a high degree of oxygenoxygen covalency.
The six valence O-O interactions in each Si(OSi) 4 MO can be categorized as bonding or anti-bonding 3 . These are illustrated and tabulated in Figure 4A  Collectively, there is a definitive excess of bonding interactions providing the net positive overlap that allows-or perhaps mandates-oxygenoxygen covalent bonding. Accordingly, the O-O bond orders are decidedly greater than zero (0.07, 0.11, 0.13, 0.13, 0.14, 0.14) and average 0.12. The Si-O bond orders are decidedly less than unity (0.88, 0.88, 0.80, 0.79) and average 0.84.
The proclamation that "there are no chemical bonds, only bonded interactions," (Gibbs et al., 2009) compels scrutiny of the Figure 4A MOs for electron density associated with O-O interactions. Figure 4B illustrates HOMO-n slices intersecting the upper right O-O axis of the O 4 tetrahedron (bond order = 0.14). These planes contain Si or are perpendicular to that-depending on which shows greater overlap. Eleven of these sixteen interactions are bonding and five are anti-bonding; locally, this O-O relationship has be O-O = 38%. Of the bonding interactions, nine utilize σ overlap and two utilize π overlap. Of the anti-bonding interactions, four are σ* and one is π*. The HOMO-13 through HOMO-9 slices reveal considerable O-O covalency, arising from three σ interactions and two π interactions. These MOs average an oxygen atomic orbital contribution of 76%-well above the 55% predicted by purely proportional silicon/oxygen atomic orbital mixing. This excess oxygen MO character portends high O-O covalency (bond order > canonical value of 0) and low Si-O covalency (bond order < canonical value of 1). The bonding excess and the abundance of O-O electron density, implied by Ψ 2 of the HOMO-n slices, are further evidence of long bonded interactions in locations not previously described.
Every atom in quartz silica is part of the Si(½O) 4 repeat unit. Additionally, every atom belongs to a twelve-membered ring of alternating silicon and oxygen atoms, Si 6 O 6 . It is illustrative to first consider the bonding arrangement of a geometry-optimized Si 6 O 6 ring, which adopts a planar, D 6h conformation. Figure 5A depicts the twenty-four valence MOs built primarily from the six oxygen 2s orbitals and the eighteen oxygen 2p orbitals (IsoValue = 0.005 [e/bohr 3 ] 0.5 ). Of the 144 O-O interactions, 68 are bonding and 76 are anti-bonding; there are 8 excess anti-bonding interactions. Hence, the bonding excess parameter is negative: be = −6%. Accordingly, the computed O-O bond orders are also negative at −0.05 (Mulliken, 1955c). Note that the HOMO-16, HOMO-13, HOMO-12, HOMO-11, HOMO-10, and HOMO-7 mimic the six π MOs of benzene (also D 6h ), built from six C 2p z atomic orbitals (Hoffmann, 1963). They possess the same number of π nodal planes (0, 1, 1, 2, 2, 3), mandating that half are bonding and half are anti-bonding. Naturally, there is no net π bonding when these six MOs are filled.
However, the Si 6 O 6 ring of quartz silica is not planar. Instead, it contorts to a chiral, C 2 -symmetric saddle shape, yielding a dramatically different array of oxygen-based valence molecular orbitals, as illustrated in Figure Figure 6B). Inspection of these silica-based molecular orbitals reveals bonding, node-free, contiguous π MOs built mostly from oxygen atomic orbitals (78%-91%, averaging 84%) and minimally from silicon atomic orbitals (9%-22%, averaging 16%). Arguably, this constitutes a modified case of homoconjugation, which is defined as π overlap separated by a single non-conjugating atom or group (McNaught and Wilkinson, 1997). In the present case, the annular π-system spans multiple non-conjugating silicon atoms, being built from atomic orbitals of relatively distant (2.61-2.64 Å) and alternating oxygen atoms. Such an arrangement that nonetheless yields bonding molecular orbitals is herein defined as alternoconjugation (Hoffmann et al., 1970).
The HOMO-17 of Figure 6B is a twist-free π MO composed mostly (83%) of six O 2p y orbitals pointing towards the middle of the Si 6 O 6 silica ring. However, the HOMO-13 and HOMO-12 are contiguous π MOs with two half-twists and thus, have a topological linking number of L k = 2 (Rappaport and Rzepa, 2008). Moreover, these MOs twist in opposite directions, resulting in quasi-enantiomorphic MOs (but not strictly enantiomorphic because the Si 6 O 6 ring itself is chiral). The HOMO-11 and HOMO-8 are also contiguous π MOs, but with four half-twists and thus, a linking number of L k = 4. These too are quasi-enantiomorphic because they twist in opposite directions. The HOMO-6 is apparently not a contiguous π MO because of a nodal surface bisecting the ring. Molecular orbitals with L k = 2 or 4 are Möbius aromatic, although such twisting generally results in a smaller resonance energy (stabilization) compared to nontwisted (Hückel) analogues with L k = 0 (Rzepa, 2005). For the HOMO-13, this Möbius aromaticity is further visualized by the six contiguous O 2p y -O 2p y π bonding interactions of Figure 6C; these HOMO-13 slices are oblique but nonetheless lack a disruptive nodal plane. Figure 6 focuses on O 2p y atomic orbitals, but there are also MOs built from O 2p z atomic orbitals ( Figure 5B and Figure 7B) that are Möbius aromatic: HOMO-15 and HOMO-14 with L k = 2; HOMO-10 and HOMO-9 with L k = 4. In total, eight of the 24 oxygen-based valence MOs are Möbius aromatic (HOMO-15 through HOMO-8). The HOMO-17 (O 2p y ) and HOMO-16 (O 2p z ) are twist-free π MOs with L k = 0 and are thus Hückel aromatic.
Altogether, the Si 6 O 6 ring of silica has ten total oxygen-based, contiguous π MOs ( Figure 5B, HOMO-17 through HOMO-8), while the D 6h Si 6 O 6 ring has only two ( Figure 5A, HOMO-17 and

Molecular orbital analysis of silica vs. sulfur
α-Sulfur is a stable chalcogen allotrope that consists of S 8 rings with sulfur-sulfur bonding (Rettig and Trotter, 1987). Compared to the O-O bonds of silica (average = 2.63 Å), the S-S bonds are shorter (average = 2.05 Å) and sulfur shares no electrons with an electropositive element. Hence, the S-S Mulliken bond orders are typical of canonical single bonds, ranging from 1.13 to 1.18. A valence molecular orbital analysis of cyclic S 8 ( Figure 7A) reveals that its S-S bonds exist because of p atomic orbital overlap and not s atomic orbital overlap-akin to the O-O bonds of the silica Si 6 O 6 ring. For the eight S 3s-based valence MOs, there are 32 bonding and 32 anti-bonding S-S interactions (be 3s = 0%). For the sixteen S 3p-based valence MOs, there are 128 bonding and 0 anti-bonding S-S interactions (be 3p = 100%). Collectively, be = 67%, allowing for net S-S bonding. Furthermore, as shown in Figure 7A and Supplementary Figure  S5, ten of the sixteen 3p-based S 8 MOs are aromatic, including: Hückel aromaticity for HOMO-12 (3p y ) and HOMO-7 (3p z ), with L k = 0; Möbius aromaticity for HOMO-9 (3p y ), HOMO-8 (3p y ), HOMO-6 (3p z ), HOMO-5 (3p z ), HOMO-2 (3p z ), and HOMO (3p z ), with L k = 2; and Möbius aromaticity for HOMO-4 (3p y ) and HOMO-3 (3p y ), with L k = 4. Of the   Figure 7B). The chalcogen-chalcogen bonding in S 8 has never been disputed. The chalcogen-chalcogen bonding in silica's Si 6 O 6 ring is equally evident and valid, given the molecular orbital and bonding parallels between S 8 and the O 6 ring embedded within silica.

Electron density analysis
The "longest O-O bond in any known molecule" was reported for HOON, studied in a supersonic molecular beam. Its O-O distance was measured to be 1.91 Å (Jabłoński, 2012;Crabtree et al., 2013;Chen et al., 2016). A computational analysis ( Figure 8A Bonding can also be gauged by the valence electron density (e/bohr 3 ) along each O-O axis, as plotted in Figure 9A. Excluding core electron density 5 the average valence electron density computes to VED ave = 0.46, 0.44, and 0.35 e/bohr 3 for HOON, O 3 , and the Si21 cluster, respectively. The valence electron density can also be parameterized according to the area under the curves of Figure

Mulliken bond orders vs subsequent methods
Although Mulliken's overlap population method "permits one to identify chemically bonded atoms," it does not provide integral bond multiplicity values (corresponding to single, double, or triple bonds) (Mayer, 2003) and it is accompanied by other various objections (Reed et al., 1985). Historically, bond orders have been calculated via several generalized methods, including those described by Mulliken (1955b), Weinhold and Landis (2005), Wiberg (1968), Löwdin (1970), Mayer (1983), and Bader (1990)-although "none of them is the "right" one." (Lewars, 2008) A comprehensive comparison is beyond the scope of this manuscript. However, a preliminary analysis demonstrates that the Mulliken, Wiberg, and Mayer methods yield comparably large O-O bond orders within silica clusters. These computed bond orders vary by only ±0.01 for a given silica cluster (see Supplementary Table  S32) (Gorelsky and Lever, 2001) 7 . Furthermore, multicenter oxygenbased bond order indices have been computed and are proportional to silica cluster size. For the large Si11-Si35 silica clusters, the 3-center 2-elecron (3c2e) bond order index ranges from I OOO = 0.319 to 0.383; these values exceed that of I BHB = 0.253 computed for diborane (B 2 H 6 ) as well as the theoretical 7 and computed value (Bochicchio et al., 1998) of 8/27 ≈ 0.296 for H 3 + . For the same large clusters, the 4-center 2electron (4c2e) bond order index ranges from I OOOO = 0.110 to 0.209; although such 4-center bonds are rare (de Giambiagi et al., 1997;Ponec and Mayer, 1997), these values mostly exceed that for B 4 (I BBBB = 0.118) (Sannigrahi and Kar, 1999) and the value of I OSiO values greater than 0.02 and only three clusters have 4c2e I OSiOO values greater than 0.01, ranging from 0.023 to 0.141. This multicenter bonding analysis highlights the electron delocalization among oxygen atoms at the expense of silicon atom participation. An additional informative metric is the atomic valence (AV). For silicic acid (Si1), the Mulliken and Mayer atomic valencies are essentially canonical with Si AV ≈ 4.33 and O AV ≈ 2.05. However, in accord with the valence electron apportionment reported in Figure 3C, the Mulliken and Mayer atomic valencies invert for large silica clusters reaching, for example, Si AV ≈ 3.14 and O AV ≈ 4.35 for Si29-confirming the high atomic valence of oxygen markedly exceeding the canonical value of 2. The Löwdin method routinely yields smaller O-O bond orders near 0.10, invariant to the silica cluster size or the location within the cluster. The Natural Bond Order (NBO) method fails with all silica clusters because of their electronic delocalization (Goodman and Sauers, 2007). The NBO method also cannot determine long bond orders in the following molecules: the terminal oxygen atoms of ozone (Hay et al., 1975); the transannular sp 2 hybridized carbons of norbornadiene (Hoffmann et al., 1970); and the transannular sulfur atoms of cyclic S 4 N 4 (Bridgeman et al., 2001). For silica, a substantial Mulliken O-O bond order is required to accord with the resonance scheme of Figure 2B (BO = 1/3). Without such long covalent bonds, Earth's most abundant material fails to have a suitable resonance formulation. Not only did Mulliken deem his method "a good measure of covalent bonding," (Mulliken, 1955b) he championed its "obvious" advantages for the interrogation of resonance and delocalized structures (Mulliken, 1955a).

Mulliken bond orders vs. computational methods and basis sets
Mulliken populations are often "unduly sensitive to basis set." (Reed et al., 1985) Table 1 reports the Mulliken bond order sensitivity for the Si5, Si11, and Si21 silica clusters subjected to various computational methods and basis sets. From the center of each cluster, the four Si-O bond orders and the six O-O bond orders are averaged to provide the values of Table 1. For the 6-311++G** basis set (used throughout this study), the bond order results are largely invariant to the computational method. Si-O bond orders are substantially less than unity and O-O bond orders are substantially greater than zero-converging as the cluster size increases and mimicking the plot for B3LYP/6-311++G** in Figure 3B.  (Mulliken, 1955c). Without diffuse functionals in the basis set (see Supplementary Table S30), the Si-O bond order is near unity and the O-O bond order is near zero, yielding the simple Lewis structure having no option for resonance. Frontiers in Chemistry frontiersin.org

Computed oxygen-oxygen bond dissociation energy for silica
If oxygen-oxygen bonding in silica is real, then perhaps an apt computational analysis could reveal silica's O-O bond dissociation energy (BDE). Conventionally, the reaction of Figure 10A between Si(OH) 4 and 3 SiH 4 is an isodesmic (Ponomarev and Takhistov, 1997) and homodesmotic (Wheeler, 2012) Figure 10D illustrates the conventionally isodesmic deconstruction of the Si29 cluster. The dendriform array of 174 O-O bonds stabilizing this cluster engenders a new term to describe the pervasive O-O bonding that supplements the conventional Si-O bonding: dendrivalent. The 8 bonds emanating from a single oxygen atom (six to oxygen and two to silicon) simultaneously depict the bonds of multiple resonance forms, none being more than trivalent at oxygen ( Figure 2B). This number of bonds will not alarm anyone fluent in resonance formulations. Such multivalent bonding is commonplace, as in the 6 bonds emanating from carbon in the resonance hybrids of the carbonate dianion (3 sigma and 3 pi) or the 9 bonds emanating from carbon in the resonance hybrids of polycyclic aromatic hydrocarbons or graphite (3 sigma, 3 pi, and arguably 3 transannular bonds of the Dewar benzene type (Sorella et al., 2014;Pauling and Wheland, 1933); 8 see Supplementary Figures S19-S21 for these multiple resonance hybrids). Inspection of valence electron density surfaces for Si5 or Si21 reveals clear valence bond paths between oxygen and all six of its oxygen neighbors ( Figure 9B; Supplementary Figures S14-S17). As noted above, the outer spherical node of the Si 3s orbital precludes silicon-oxygen valence bond paths (Figure 9, VED min = 0.000 e/bohr 3 ).
While silica employs more O-O valence electrons than Si-O valence electrons according to Figure Figure 10F. The prediction is [114.2 + 4.4 + 4.4 + 4.4] kcal/mol or 127.4 kcal/mol, which is very close to the enthalpy computed for the homolytic bond cleavage of HO-Si(OH) 3 , ΔH = +127.6 kcal/mol.
The similarity of these values indicates that a consistent and predictive bond energy additivity scheme can be formulated by proper inclusion of both long and short covalent bonds. For example, the greater C-F bond strength in CF 4 (130.5 kcal/mol) versus that in CH 3 F (109.9 kcal/mol) has been rationalized by electronegativity and atomic charge effects (Lemal, 2004). However, a bond energy additivity scheme, analogous to that for silica in Figure 10, provides a simpler explanation. The computational isodesmic reaction of CF 4 + 3 CH 4 → 4 CH 3 F has ΔH = +39.4 kcal/mol ( Figure 10G), suggesting an F-F BDE of 6.6 kcal/ mol since this reaction breaks six long (2.17 Å), unconventional F-F bonds. Separately, the manifold C-F BDE in F-CF 3 , which represents one C-F bond and three F-F bonds, is predicted by adding the experimental solitary C-F BDE in F-CH 3 to three F-F BDE values: [109.9 + 6.6 + 6.6 + 6.6] kcal/mol = 129.7 kcal/mol. This result is strikingly close to the experimental value of 130.5 kcal/mol and does not invoke nebulous arguments about electronegativity/ atomic charge (Lemal, 2004), negative hyperconjugation (Hine, 1963;Reed and von Ragué Schleyer, 1987), or Coulombic interactions (Wiberg and Rablen, 1993), but elaborates on double bond/nobond resonance formulations (Dolbier et al., 1982) by adding six new resonance hybrids, each with a long F-F bond, as in [F-C-F][F-F] (explicitly drawn in Supplementary Figure S22). A valence molecular orbital analysis corroborates the F-F bond in this unconventional resonance hybrid for CF 4 since be 2s = 0% and be 2p = 39%, yielding a composite bonding excess greater than zero at be F-F = 29%, which is rather close to be O-O = 27% for the silicic acid model Si(OSi) 4 (Figure 4), or be O-O = 21% for silicic acid itself. To recapitulate: CF 4 and silica clusters have a substantial excess of bonding versus antibonding F-F or O-O interactions that mandate an F-F or O-O bond dissociation energy of non-zero magnitude. While these long BDE values are small compared to those of short bonds, there are six per tetrahedron and thus, they stabilize CF 4 by 39.4 kcal/mol and silica

Ball-and-stick model vs cooperative sphere model
For silica, it is evident that silicon-oxygen bond dissociation energies are cooperative (Figures 10E, F) and there are no valence bond paths between oxygen and silicon because of an intervening nodal surface ( Figures 9A, B). Hence, a ball-and-stick bonding model-germane to water, methane, and many other molecules-is too simplistic to accurately describe the bonding in silica. For Si[OSi(OH) 3 ] 4 (the Si5 silica cluster), the oxygen and silicon contributions to the valence electron density can be functionally separated into the 2D plots of Figure 11A 9 .
Electron density from the silicon 3s and 3p atomic orbitals accumulates about 0.8 Å from silicon, coinciding with the midpoint of the geminal O-O axes, where oxygen 2s and 2p atomic orbitals also accumulate electron density. Thus, silicon and oxygen are perfectly matched to create a sphere of electron density 0.8 Å from silicon, about halfway along the silicon-oxygen axes, which are 1.61 Å long; Figure 11A locates this sphere (in orange) just beyond the outer spherical node of the Si 3s atomic orbital, which encircles silicon 0.4 Å away (where VED = 0 in Figure 9A). Because of this nodal surface, the silicon-oxygen bond of silica does not merit a canonical line between the nuclei; instead, bonding is better defined by a cooperative sphere model, where the valence electron density accumulates in a spherical region 0.8 Å from silicon and 0.8-1.3 Å from oxygen.
Between each pairwise O-O interaction of the Si5 central silicate, the valence electron density bond paths reach minima of 0.026, 0.026, 0.026, 0.026, 0.025, and 0.025 e/bohr 3 ( Figure 11A). Notably, the same analysis computed without diffuse functionals ( Figure 11B) yields the exact same minima between each pairwise O-O interaction and generates indistinguishable valence electron density 3D surfaces and 2D plots. Both diffuse and non-diffuse computational methods aspire to the same valence electron density distribution-with equivalent accumulation of electron density between all geminal oxygen atoms and a nodal surface between canonically bonded silicon and oxygen atoms. However, the two computational methods employ oxygen and silicon atomic orbitals differently to

FIGURE 11
For the Si5 silica cluster, the computed valence electron density distribution with diffuse functionals (A) is nearly identical to that without diffuse functionals (B), despite differential employment of the constituent oxygen and silicon atomic orbitals. This differential accounting of each electron's orbital origin results in non-canonical (A) or canonical (B) O-O and Si-O bond orders. Nonetheless, both computational methods show a spherical node between silicon and oxygen, but a definitive accumulation of electron density between all geminal oxygens, resulting in a sphere of electron density (orange circle) responsible for the cooperative bonding in silica.
9 The valence electron density of Si5 was parsed to the oxygen valence electron density component by setting all silicon and hydrogen coefficients to zero in the Spartan Output (.txt) file prior to the Analyze MO compositions function in Chemissian. The silicon valence electron density component was similarly computed after setting all oxygen and hydrogen coefficients to zero.
Frontiers in Chemistry frontiersin.org build the valence molecular orbitals. With diffuse functionals, the valence molecular orbitals are 70% oxygen-based and 27% siliconbased. This allocation is more balanced than the computation without diffuse functionals, for which the valence molecular orbitals are 81% oxygen-based and 15% silicon-based. (In both cases, about 3% of the valence molecular orbitals are hydrogen-based.) A purely canonical model intervenes with valence molecular orbitals that are 75% oxygen-based (96 e/128 e) and 16% silicon-based (20 e/128 e). Furthermore, with diffuse functionals, the oxygen component shows significant overlap of the oxygen atomic orbitals ( Figure 11A) Since the diffuse and confined computational models allocate electrons differently to atomic orbitals, they differ in the resultant bond orders. Nonetheless, the valence electron density is unmistakably similar, as exhibited by the 3D and 2D valence electron density plots of Figure 11 as well as the 1D plots of Figure 12 for the Si21 silica cluster. Specifically, the VED, VED ave , and VED proj along O-O or O-Si axes are essentially invariant to the employment of diffuse functionals. But, regardless of whether the computational model includes or excludes diffuse functionals, more valence electrons reside between oxygen atoms (VED ave = 0.35; VED proj = 1.73 or 1.72) than between silicon and oxygen atoms (VED ave = 0.33; VED proj = 0.99).

Overlap population density of states analysis
Additional bonding insight for silica is provided by the valence density of states map (Lu and Chen, 2012) shown in Figure 13. For a central oxygen-oxygen interaction in the silica clusters Si1, Si5, Si11, Si15, Si21, and Si25, the partial density of states (PDOS, orange thin lines) is plotted as a function of molecular orbital energy. Also plotted is the overlap population density of states (OPDOS) (Hughbanks and Hoffmann, 1983), which reveals the bonding states (positive) and anti-bonding states (negative) for a series of molecular orbitals-in analogy to the crystal orbital overlap population (COOP) (Grechnev et al., 2003)  In the molecular orbital region of −29 to −26 eV for all six clusters, the OPDOS curves ( Figure 13, thick blue lines) are both positive and negative. This suggests that MOs built from O 2s atomic orbitals do not provide substantial net bonding. However, in the molecular orbital region of −18 to −8 eV, the OPDOS curves become increasingly positive with increasing silica cluster size. This suggests that O 2p atomic orbital overlap is largely responsible for the oxygen-oxygen bonding in silica. Exemplary molecular orbitals are shown in Figure 13 for the Si11 cluster. Oxygen atom #1 and oxygen atom #2 have σ bonding interactions via their 2s atomic orbitals in the HOMO-125, their 2p z atomic orbitals in the HOMO-93, and their 2p x/y atomic orbitals in the HOMO-35. Orthogonal cross sections of these MOs locate the accumulation of electron density along the O-O axis. Importantly, the OPDOS analysis shows that increasing cluster size augments the bonding for the O 2p x/y -O 2p x/y interactions (−12 to −8 eV) more than the O 2p z -O 2p z interactions (−18 to −12 eV).
What is the basis for this bonding augmentation? A detailed OPDOS analysis of all six central O-O valence interactions in the Si1-Si25 clusters (Table 2) reveals two effects. The first effect concerns the number of bonding vs anti-bonding interactions. Adding oxygen atoms to the periphery of a silica cluster decreases the nodal density between oxygen atoms in the center Frontiers in Chemistry frontiersin.org of the silica cluster. The nodal density is inversely related to the OPDOS % bonding excess reported in Table 2. While the % bonding excess is somewhat invariant to cluster size for the O 2s MOs (−29 to −26 eV) and the O 2p z MOs (−18 to −12 eV), this parameter increases correspondingly with cluster size for the O 2p x/y MOs (−12 to −8 eV)-from −38% to 26%-fully consistent with the dominant growth of the OPDOS curve in the same region ( Figure 13). The second effect concerns the magnitude of the bonding vs anti-bonding interactions. Adding oxygen atoms to the periphery of a silica cluster increases the cumulative OPDOS bonding parameter more than the cumulative anti-bonding OPDOS parameter for the oxygen atoms in the center of the silica cluster. Thus, the |bonding/anti-bonding| quotient steadily increases from 0.76 to 2.42 with increasing cluster size (Table 2).
(Note that the |bonding/anti-bonding| quotient in Figure 13 represents only one of the six O-O interactions in the central silicate unit of the clusters-all of which are averaged to yield the values in Table 2

Conclusions and outlook
A molecular orbital analysis of α-quartz silica model complexes reveals that oxygen valence electrons abandon their canonically prescribed locations to form long covalent oxygen-oxygen bonds. Oblique arrangements of the oxygen atoms minimize molecular orbital nodes and maximize bonding interactions that are 2.61-2.64 Å apart, yet have Mulliken bond orders reaching 0.63 and averaging 0.47. This nodal minimization is inherently challenging for linear or planar molecules, but may prove widespread among atoms arranged in three dimensions-accomplished especially via p atomic orbitals. Thereby, the geminal oxygen atoms of silica bond cooperatively with a computed O-O bond dissociation energy of 4.4 kcal/mol Frontiers in Chemistry frontiersin.org when they are suitably arranged and interspersed with the electropositive element silicon. Generally, as silica model complexes increase in size, canonical bonding paradigms decrease in accuracy and resonance hybrids increase in relevance. This comports with the original Hund-Mulliken molecular orbital theory, which states that the "best MOs . . . spread at least to some slight extent over all atoms." (Mulliken, 1970) For Pauling, the concept of resonance had "its most important chemical applications" to "molecules to which no satisfactory single structure in terms of single bonds, double bonds, and triple bonds can be assigned." (Pauling, 1946) Indeed, the structure, bonding, and energetic stability of silica cannot be fully understood without resonance hybrids involving long oxygen-oxygen bonds. This covalent bonding likely exists between other distant atoms and promises to impact the understanding of many materials and processes. The hybridization theory of Pauling (Pauling, 1930;Pauling, 1931) and Slater (Slater, 1931) compelled chemists to reconsider the location of nearly all valence electrons. The curious bonding in ferrocene (Fischer and Pfab, 1952;Wilkinson et al., 1952), with multiple metal-carbon bonds, also provided a paradigm shift in our understanding of bonding and the location of electrons; however, this bonding arrangement pertains only to a class of esoteric, manmade organometallic sandwich compounds. The bonding in diborane (Lipscomb, 1973), with three-center two-electron bonds, provided another paradigm shift in our understanding of bonding arrangements in molecules; but this bonding motif is also uncommon, primarily ascribed to a limited number of main group and transition metal complexes (Longuet-Higgins and Roberts, 1955). The long covalent bond theory (LCBT) posited herein signals a new paradigm shift in the location of the chemical bond. As evidenced by resonance formulations, bond orders (Mulliken, Wiberg, Mayer), multicenter bond order indices (3center and 4-center), atomic valencies (Mulliken, Mayer), molecular orbitals (bonding > anti-bonding interactions), bonding analogies, electron density calculations, valence bond path calculations, overlap population density of states (OPDOS) analysis, and bonding energetics, it is clear that Nature builds materials with long covalent bonds-not just the short canonical bonds introduced by Couper (1858) and Kekulé (1866) whose primitive bonding model has somehow reigned over 160 years. These long bonds underpin an explanation for the contorted, chiral structure and energetic stability of α-quartz silica, wherein distant oxygen-oxygen bonding supplements conventional bonding. This bonding paradigm is abundant and pervasive since more than onethird of the 10 49 valence electrons in the Earth's crust are allocated to long covalent bonds 10 . Moreover, 10 48 crustal valence electrons inhabit silica's oxygen-based Möbius aromatic orbitals-manifestly the most prevalent sort of aromaticity. Astonishingly, LCBT implicates the oxygen-oxygen bond as the most abundant bond on Earth. While this study focuses on silica, future work will reveal the prevalence, energetics, and importance of long covalent bonds in a rather wide variety of materials-especially those with periodic structures, including ice, biopolymers, bone, and superconducting ceramics. 10 Silica constitutes 59% of the Earth's crust (Clarke and Washington, 1924) and provides about 66% of crustal valence electrons. Since over half (52%) of these electrons are spent on oxygen-oxygen bonding, over one-third of crustal valence electrons participate in the long O-O covalent bonds of silica.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.