Planar tetracoordinate oxygen stabilized within triel–chalcogen dicationic frameworks

Abstract

Planar tetracoordinate oxygen (ptO) challenges oxygen's intrinsic preference for two-coordinate bonding. We present a systematic investigation of dicationic triel–chalcogen clusters O©X₄E₄²⁺ (X = Al, Ga, In; E = S, Se, Te), identifying 9 D_4h ptO global minima. The oxygen atom is best described as a formally O²⁻ dianion electrostatically embedded within a charge-delocalized X₄E₄ tetracationic framework sustained by a robust X–E σ-scaffold, with negligible X–X interactions. Bonding analysis indicates predominantly electrostatic O–X interactions, whereas the X₄E₄ crown is stabilized by localized 2c–2e X–E σ-bonds and delocalized 3c–2e X–E–X π-bonds. Magnetic response calculations reveal localized diatropic vortices around the ptO center and the X–E–X regions, but no global ring current, consistent with the nonaromatic character. Born–Oppenheimer molecular dynamics simulations confirm dynamic stability at 500 K. π-Anion complexation with cyclooctatetraenide (COT²⁻) demonstrates that retention of the planar geometry depends on intrinsic framework rigidity and is preserved in the aluminum derivatives. These results establish electrostatic embedding within a rigid X₄E₄ scaffold as an effective strategy for stabilizing ptO in main-group clusters.

---

1 Introduction

The concept of planar hypercoordination has significantly broadened classical bonding paradigms in main-group chemistry, challenging conventional valence expectations and revealing new modes of electronic stabilization.¹ While early theoretical work on planar tetracoordinate carbon established that tetrahedral preferences can be overcome through appropriate electronic compensation and charge redistribution, extending this concept to oxygen remains particularly demanding.^2–10 Owing to its high electronegativity and strong tendency toward localized two-center bonding, oxygen intrinsically disfavors planar hypercoordination.

Planar tetracoordinate oxygen motifs have been experimentally identified in molecular and solid-state systems, where the oxygen atom is embedded within extended electropositive environments such as metal clusters and coordination frameworks.^11–13 In addition, theoretical studies have identified simple group 13 systems such as Al₄O as prototypical ptO motifs. In this context, the O©Al₄ unit investigated here can be viewed as an idealized molecular model that captures the essential electrostatic stabilization features of these environments.^14,15 A major advance was the systematic potential energy surface (PES) screening studies by Merino and co-workers,¹⁶ which identified 35 ptO clusters, rigorously validated as true global minima (GMs), predominantly stabilized by Group 13 elements. By establishing these species as true GMs on their PESs, these studies showed that ptO can be thermodynamically preferred rather than kinetically trapped, and delineated stabilization regimes ranging from covalent to strongly polar, charge-assisted interactions. Complementary developments further broadened the design space: Li et al. proposed noble-gas-supported ptO dications as GMs, emphasizing polarization and weak-interaction contributions in unconventional bonding environments; however, subsequent PES analyses indicate that these structures do not correspond to true GMs.^12,13 Additionally, Wang et al. proposed OLi₄Au₄²⁻ as a GM structure in which ptO is stabilized through dominant coulombic confinement within a rigid peripheral framework.¹⁷ Notably, planar hypercoordination at oxygen extends beyond tetracoordination. Wu et al. demonstrated that the gas-phase [Be₅O₆]²⁻ dianion, previously misassigned as linear, in fact possesses a planar pentacoordinate oxygen (ppO) GM, underscoring the importance of rigorous structural characterization and comprehensive PES exploration in validating planar hypercoordinate motifs.^18,19

Despite the identification of numerous ptO global minima across diverse compositions, stabilization within a symmetry-preserving and compositionally tunable framework remains insufficiently examined. We therefore focus on a D_4h X₄E₄ platform that enables systematic triel and chalcogen substitution within a fixed topology. The O©X₄E₄²⁺ framework builds on a structural motif previously employed in chalcogen-bridged clusters featuring planar hypercoordinate centers, such as N©X₄E₄⁺ systems.^20–22 In these architectures, the X₄E₄ ring is defined by a combination of localized X–E σ bonds and delocalized 3c–2e X–E–X interactions, which together provide electronic compensation of the electron-deficient metal centers and confer substantial structural rigidity. This bonding pattern stabilizes electronegative central atoms in a planar environment. In this context, we extend this design principle to oxygen by embedding a formally O²⁻ center within a preorganized D_4h–X₄E₄⁴⁺ scaffold, preserving the characteristic bonding topology while enabling stabilization of a planar tetracoordinate oxygen center (Scheme 1). This construction yields a closed-shell 40-valence-electron architecture in which ptO stabilization arises from charge compensation and geometric confinement.

Scheme 1 Electrostatic embedding of a formally O²⁻ center within a preorganized D_4h–X₄E₄⁴⁺ scaffold yields a closed-shell O©X₄E₄²⁺ complex with preserved topology.

To evaluate this design hypothesis, we systematically examined the dicationic O©X₄E₄²⁺ (X = B–Tl; E = O–Po) clusters, featuring a central oxygen embedded within a square triel array bridged by chalcogen ligands. Among 25 symmetry-constrained D_4h candidates, geometry optimizations and vibrational analyses locate 16 stationary ptO minima. Exhaustive global PES exploration establishes that 9 correspond to true GM on their respective PES, validating this triel–chalcogen framework as a tunable and symmetry-preserving ptO platform. A representative example is D_4h-O©Al₄S₄²⁺, in which oxygen remains strictly planar within a square Al₄ unit supported by sulfur bridges.

To assess stabilization beyond the isolated gas phase, complexation with the aromatic COT²⁻ dianion was examined. While the resulting neutral assemblies correspond to minima on the PES, preservation of the planar ptO core is observed only for the aluminum derivatives. In contrast, heavier triel systems undergo partial distortion of the X₄E₄ framework due to stronger covalent interaction with the electron-rich π surface of COT²⁻. These results indicate that counterion-assisted stabilization of ptO motifs is feasible when the interaction remains predominantly electrostatic, whereas excessive covalent engagement compromises planarity.

2 Computational methods

The PESs of the O©X₄E₄²⁺ systems (X = Al, Ga, In; E = S, Se, Te, Po) were explored using the AUTOMATON program,²³ which employs a probabilistic cellular automata algorithm combined with genetic operations to generate and evolve candidate structures toward GM. Initial screening was performed at the PBE0 (ref. 24)/SDDAll^25–29 level for both singlet and triplet states. Structures lying within 90 kcal mol⁻¹ of the lowest-energy isomer were subsequently reoptimized at the PBE0 (ref. 24)-GD3 (ref. 30)/def2-TZVP³¹ level of theory. Harmonic vibrational frequency calculations were carried out to confirm that all retained structures correspond to true minima on their respective PESs. For each identified GM, wavefunction stability was verified through single-point calculations using the stable = opt procedure, confirming that the self-consistent field (SCF) solution corresponds to a minimum in wavefunction space. Final energetic refinements were obtained from single-point DLPNO-CCSD(T)^32–34/CBS^35,36//PBE0-GD3/def2-TZVP calculations, including zero-point energy (ZPE) corrections. All geometry optimizations, frequency calculations, and wavefunction stability analyses were performed using Gaussian 16.C02,³⁷ whereas DLPNO-CCSD(T)/CBS single-point calculations were carried out with ORCA 6.0.1.³⁸ To assess the single-reference character of the electronic structures, the T₁-diagnostic³⁹ was computed from the DLPNO-CCSD(T) wavefunctions.

Bonding analyses were conducted using Natural Population Analysis⁴⁰ (NPA) and Wiberg bond indices⁴¹ (WBI) within the Natural Bond Orbital (NBO 7.0) framework.^42,43 The Adaptive Natural Density Partitioning (AdNDP) method,⁴⁴ as implemented in Multiwfn 3.8,⁴⁵ was employed to identify both localized and multicenter bonding elements. Interatomic interactions were further characterized using the Interacting Quantum Atoms^46–48 (IQA) formalism at the PBE0-GD3/def2-TZVP level with the AIMAll program,^49–52 where the interaction energy (V^int_IQA) is expressed as the sum of coulombic (V^int_C) and exchange–correlation (V^int_XC) contributions.

The possible aromatic character of the D_4h–O©X₄E₄²⁺ clusters was examined through magnetically induced current density⁵³ (MICD) calculations using the GIMIC program^54,55 at the ωB97X-D⁵⁶/def2-TZVP³¹ level. The use of a range-separated hybrid functional is motivated by previous studies showing that long-range Hartree–Fock exchange is essential for a realistic description of electron delocalization in extended π-conjugated and macrocyclic systems, where conventional hybrid functionals may overestimate cyclic delocalization.⁵⁷ An external magnetic field was applied perpendicular to the molecular plane (z-axis), and the gauge-including atomic orbital⁵⁸ (GIAO) formalism was employed to ensure origin-independent results. To assess the extent of delocalization in the investigated systems, Electron Density of Delocalized Bonds (EDDB)^59,60 calculations were performed at the ωB97X-D/def2-TZVP level of theory. Previous studies have shown that EDDB results are insensitive to basis set size.⁶⁰

3 Results and discussion

3.1 Systematic PES screening and identification of minimum-energy ptO structures

As an initial screening step, the 25 O©X₄E₄²⁺ combinations (X = B, Al, Ga, In, Tl; E = O, S, Se, Te, Po) were fully optimized under D_4h symmetry at the PBE0²⁴ -GD3 ³⁰/def2-TZVP³¹ level of theory, followed by harmonic frequency analysis to confirm their nature as stationary points. 16 systems were verified as genuine planar minima on the PES, whereas the remaining cases either relaxed to non-planar structures or corresponded to higher-order saddle points (Scheme 2 and Table S1). Notably, imaginary frequencies are restricted to the boron derivatives and to O©X₄O₄²⁺series, whereas all Al, Ga, and Tl chalcogen systems exhibit positive lowest vibrational frequencies, confirming their intrinsic planar stability at the DFT level.

Scheme 2 Workflow allowing the identification of global minima of O©X₄E₄²⁺.

These 16 planar candidates were subsequently subjected to an unbiased global exploration of their PES using the AUTOMATON algorithm²³ (see Computational Details). This second-stage search revealed 9 ptO GMs on their respective PES, namely O©Al₄S₄²⁺, O©Ga₄S₄²⁺, O©In₄S₄²⁺, O©Al₄Se₄²⁺, O©Ga₄Se₄²⁺, O©In₄Se₄²⁺, O©Al₄Te₄²⁺, O©Ga₄Te₄²⁺, O©In₄Te₄²⁺. For all nine confirmed ptO GMs, the D_4h structure (labelled “a” in Fig. S1–S9) is clearly separated from the second-lowest isomer. At the DLPNO-CCSD(T)^32–34/CBS^35,36//PBE0-GD3/def2-TZVP level of theory (including ZPE corrections), the energy differences between the global minimum and the second isomer (ΔE) amount to 61.6 kcal mol⁻¹ for O©Al₄S₄²⁺, 50.7 kcal mol⁻¹ for O©Ga₄S₄²⁺, 5.6 kcal mol⁻¹ for O©In₄S₄²⁺, 60.4 kcal mol⁻¹ for O©Al₄Se₄²⁺, 60.4 kcal mol⁻¹ for O©Ga₄Se₄²⁺, 18.0 kcal mol⁻¹ for O©In₄Se₄²⁺, 50.7 kcal mol⁻¹ for O©Al₄Te₄²⁺, 50.7 kcal mol⁻¹ for O©Ga₄Te₄²⁺, and 50.7 kcal mol⁻¹ for O©In₄Te₄²⁺. Except for the indium derivatives—particularly O©In₄S₄²⁺, where the gap is modest—the ptO structures are strongly stabilized relative to competing minima, confirming that ptO corresponds to a well-defined and energetically isolated GM.

A clear compositional pattern emerges from the global exploration. Genuine D_4h ptO GMs are exclusively obtained for X = Al, Ga, and In combined with E = S, Se, and Te, whereas boron, thallium, oxygen, and polonium derivatives are not GMs. This trend reveals a narrow stabilization window requiring intermediate triel electropositivity together with a chalcogen capable of sustaining a robust yet not excessively polar X–E σ/π framework. Boron derivatives exhibit increased covalent character that weakens the electrostatic confinement of O²⁻, while thallium systems suffer from reduced electrostatic stabilization due to orbital diffuseness. Similarly, E = O will over-polarize the framework, whereas E = Po provides insufficient σ reinforcement. Thus, ptO stabilization in this platform follows a “Goldilocks” regime in which electrostatic embedding and peripheral covalent reinforcement are optimally balanced.

The electronic structures of these 9 GMs further support their robustness (Table S2). The lowest vibrational frequencies of the nine global minima range from 15.3 to 67.8 cm⁻¹ (Table S1). In addition, the calculated IR spectra of all nine GMs were simulated (Fig. S10). The most intense absorption band appears in the 480–650 cm⁻¹ region and corresponds to a doubly degenerate E_u mode, fully consistent with the D_4h symmetry of the clusters. Analysis of the displacement vectors shows that this vibration does not correspond to a localized bond-stretching motion, but rather to a collective skeletal distortion involving the peripheral framework. The strong IR activity of this mode reflects the substantial dipole-moment variation generated by this cooperative motion, making it the main spectroscopic fingerprint of the O©X₄E₄²⁺ series. All systems display large HOMO–LUMO gaps (3.8–5.5 eV) and substantial singlet–triplet separations (57.7–85.7 kcal mol⁻¹), indicative of electronically closed-shell singlet ground states. Moreover, the T₁-diagnostics (≈0.01 in all cases) confirm the single-reference character of the wavefunctions, validating the reliability of the DLPNO-CCSD(T) treatment. Collectively, vibrational stability, energetic isolation, sizable electronic gaps, and negligible multireference character consistently establish these species as genuine, electronically well-defined ptO GMs.

Born–Oppenheimer molecular dynamics simulations at 500 K were performed to assess the dynamical stability of the nine D_4h GMs. Representative RMSD profiles for the aluminum series O©Al₄E₄²⁺ (E = S, Se, Te) are shown in Fig. 1, while the corresponding Ga and In derivatives are provided in the SI (Fig. S11 and S12). For O©Al₄S₄²⁺, the RMSD remains tightly bound within ∼0.08–0.20 Å over the entire 10 ps trajectory. O©Al₄Se₄²⁺ exhibits slightly larger oscillations, typically ∼0.12–0.30 Å with occasional peaks approaching ∼0.35 Å. In contrast, O©Al₄Te₄²⁺ displays pronounced initial deviations during the early equilibration stage, reaching ∼0.9–1.0 Å, followed by stabilization within ∼0.30–0.45 Å for the remainder of the simulation. The Ga and In analogues (see SI) show analogous behavior: the S derivatives remain within ∼0.10–0.20 Å, the Se systems fluctuate mainly within ∼0.15–0.30 Å, and the Te derivatives display broader but still bounded oscillations generally below ∼0.40–0.45 Å after equilibration. Importantly, in all nine trajectories, no bond rupture, isomerization, or structural rearrangement is observed, and the ptO motif is preserved throughout, confirming that the D_4h structures correspond to dynamically stable configurations under thermal perturbation.

Fig. 1 RMSD curves along the trajectories of Born–Oppenheimer molecular dynamics for a O©Al₄E₄²⁺ systems (E = S, Se, Te). The X-axis shows time in femtoseconds (fs), while the Y-axis shows RMSD in angstroms (Å). Systems E: (blue = S, red = Se and green = Te). Temperature (500 K).

3.2 Structure and bonding of the D_4h ptO global minima

The structural and electronic parameters of the 9 D_4h-O©X₄E₄²⁺ (X = Al, Ga, In; E = S, Se, Te) GMs are summarized in Fig. 2. In all cases, the central oxygen bears a substantial negative charge (q(O) = −1.49 to −1.56 |e|), whereas the group 13 atoms are strongly electropositive (q(X) = +0.95 to +1.58 |e|), confirming pronounced charge separation and significant ionic character in the O–X interactions. The peripheral chalcogens carry moderate negative charges (q(E) = −0.70 to −0.08 |e|), decreasing systematically from S to Te, consistent with reduced electronegativity down the group. The charge distribution indicates an electron-rich oxygen center stabilized within an electropositive framework, dominated by electrostatic interactions.

Fig. 2 (a) Schematic representation of the global minima detected in this work. (b) Electrostatic potential and NPA charges (q, |e|) of the O©X₄E₄²⁺ systems (c) bond lengths (in Å) and Wiberg bond indices (WBI) of the O©X₄E₄²⁺ computed at the PBE0-GD3/def2-TZVP level.

The O–X Wiberg bond indices are uniformly low (WBI_O–X = 0.17–0.20), indicating weak covalent contributions, suggesting dominant electrostatic stabilization. In contrast, the X–E bonds display substantial bond orders (WBI_X–E = 0.88–1.13), consistent with classical σ-bonding within the peripheral framework. The X–X interactions are negligible (WBI_X–X = 0.04–0.11), confirming the absence of direct metal–metal bonding and supporting a star-like, oxygen-centered bonding topology. Geometrically, r_O–X increases monotonically down group 13 (Al < Ga < In), ranging from 1.93–1.97 Å (Al) to 2.21–2.24 Å (In), mirroring the increase in atomic radius. These distances are slightly longer than the sums of Pyykkö single-bond radii (1.89, 1.87, 2.05 Å for O–Al, O–Ga, and O–In, respectively) yet significantly shorter than the corresponding van der Waals separations (3.36–3.52 Å), confirming genuine bonding interactions rather than mere electrostatic contacts. Similarly, the X–E distances increase from S to Te and remain close to the sums of single-bond radii while far below van der Waals limits, consistent with strong peripheral σ-framework stabilization.

Collectively, the data establish a bonding picture dominated by electrostatic O²⁻⋯X⁺ interactions reinforced by a robust X–E σ-framework. The planar tetracoordinate oxygen is therefore stabilized not by classical covalent hypercoordination, but by charge-assisted bonding within an electronically compensated, highly symmetric D_4h environment.

To further elucidate the electronic origin of the planar tetracoordinate oxygen motif, an AdNDP analysis was carried out using O©Al₄S₄²⁺ as a representative example in Fig. 3 (the remaining systems are depicted in Fig. S13–S15). The peripheral Al–S framework recovers four localized S lone pairs (ON ≈ 1.95–1.96 |e|), eight 2c–2e Al–S σ bonds (ON ≈ 1.94–1.97 |e|), and four 3c–2e Al–S–Al π bonds (ON ≈ 2.00 |e|), confirming a robust σ/π-stabilized crown fully consistent with the large WBI_X–E values discussed above, and with our proposed design strategy (Scheme 1). The oxygen-centered interaction admits two chemically consistent partitions. In partition scheme b (Fig. 3b), the central atom is described by four 1c–2e lone pairs with occupation numbers in the range 1.82–1.90 |e|; their deviation from the ideal 2.00 |e| reflects partial delocalization of oxygen density toward the Al₄ fragment. This representation aligns with the charge analysis (q(O) = −1.54 |e|) and the small WBI_O–Al value (0.19), supporting an O²⁻-like center electrostatically embedded within a cationic Al₄S₄ framework. In contrast, partition scheme c (Fig. 3c) describes the same region through three 5c–2e O–Al₄ σ bonds and one 5c–2e π bond, for which the occupation numbers are consistently closer to 2.00 |e| (σ: 1.96–1.99 |e|; π: 1.96 |e|), indicating a more internally coherent multicenter representation. The same qualitative behavior is observed for the Ga and In analogues and for the Se and Te derivatives (Fig. S13–S15), where q_O remains large (−1.49 to −1.56 |e|), WBI_O–X remains small (0.17–0.20), and the delocalized description systematically yields occupation numbers nearer to 2.00 |e| than the localized lone-pair scheme. Taken together, the AdNDP results support the electrostatic picture established from structural and charge analysis: the ptO center retains pronounced anionic character, yet its stabilization in the D_4h environment arises from a synergy between dominant O²⁻⋯X⁺ electrostatic interactions and secondary multicenter O–X delocalization within a strongly covalent X–E scaffold.

Fig. 3 AdNDP bonding pattern of O©X₄E₄²⁺ (X = Al–In; E = S–Te). (a) AdNDP scheme with four 1c–2e O lone pairs, eight 2c–2e E–X σ-bonds and four 3c–2e X–E–X π-bonds. (b) Alternative pattern in line with an electrostatic model. (c) Alternative pattern in line with delocalized bonding model.

To quantitatively assess the bonding picture, an IQA energy decomposition analysis was performed for the 9 D_4h GMs. Fig. 4 reports the sulfide series, which is discussed here as representative, while the corresponding data for the Se and Te analogues are compiled in Table S3.

Fig. 4 Energy components of IQA for the D_4h O©X₄S₄²⁺ systems (X = Al–In); V^int_IQA, V^int_C, and V^int_XC are interatomic IQA interaction energy and their coulombic and exchange-correlation energy components, respectively, in kcal mol⁻¹.

For the sulfides O©X₄S₄²⁺, the O–X interaction is strongly stabilizing and predominantly coulombic. The total interatomic interaction energies V^int_IQA (O–X) are −744.2 (Al), −426.2 (Ga), and −348.2 kcal mol⁻¹ (In), with the Coulomb component V^int_C (O–X) accounting for −714.9, −368.4, and −291.7 kcal mol⁻¹, respectively, whereas the exchange–correlation term V^int_XC (O–X) is comparatively small (−29.2 to −57.8 kcal mol⁻¹). Thus, across the series, O–X stabilization is mainly electrostatic, consistent with the large negative q(O) values and the low WBI_(O–X) indices (0.17–0.20). The non-negligible V^int_XC contribution nonetheless confirms a secondary covalent component, consistent with the weak multicenter O–X delocalization identified by AdNDP. X–X interactions are strongly destabilizing and essentially purely coulombic. For the sulfides, V^int_IQA (X–X) equals +661.2 (Al), +280.8 (Ga), and +215.8 kcal mol⁻¹ (In), with V^int_C (X–X) accounting for nearly the entire repulsion and V^int_XC (X–X) remaining negligible. These values quantitatively corroborate the negligible WBI (X–X) values and confirm that metal–metal bonding does not contribute to stabilization; instead, X–X electrostatic repulsion represents an intrinsic energetic penalty of the square-planar geometry.

In contrast, X–E interactions are strongly attractive and exhibit substantial exchange–correlation contributions. For O©X₄S₄²⁺, V^int_IQA (X–X) amounts to −665.5 (Al), −294.1 (Ga), and −239.9 kcal mol⁻¹, with sizable V^int_XC (X–X) terms (−67.4 to −123.8 kcal mol⁻¹), confirming covalent character. This is fully consistent with the high WBI (X–E) values and the σ/π bonding framework recovered by AdNDP. The remaining six global minima, O©X₄Se₄²⁺ and O©X₄Te₄²⁺ (see Table S3, SI), display the same hierarchy of interactions: O–X interactions remain predominantly coulombic and stabilizing, X–X interactions remain repulsive and nearly purely electrostatic, and X–E interactions retain significant exchange–correlation contributions indicative of strong covalent σ-framework stabilization. Although the absolute magnitudes decrease systematically from Al to In and from S to Te—consistent with increasing bond lengths and reduced charge density overlap—the relative partitioning of stabilization is preserved across the full 3 × 3 series.

The IQA delocalization indices (δ) further support the bonding picture. For the sulfide series, δ (O–X) range from 0.18 (Al) to 0.13 (Ga) and 0.12 (In), confirming limited electron sharing between oxygen and the group 13 centers, in agreement with the low WBI_(O–X) values (0.17–0.20). In contrast, δ (X–S) values are significantly larger, 0.78 (Al), 0.63 (Ga), and 0.58 (In), consistent with the substantial exchange–correlation stabilization and the high WBI_(X–E) indices, indicating pronounced covalent character within the peripheral framework. The δ (X–X) values remain very small (0.06–0.03), corroborating the absence of direct metal–metal bonding. The Se and Te analogues (Table S3, SI) display the same hierarchy, with δ (O–X) remaining low (≈0.10–0.17), δ (X–E) substantially larger (≈0.55–0.80), and δ (X–X) negligible. Thus, the delocalization indices quantitatively confirm that electron sharing is concentrated in the X–E scaffold, whereas the O–X interaction is primarily electrostatic with only secondary covalent character.

Collectively, the IQA results confirm and quantify the bonding model: stabilization of the planar tetracoordinate oxygen motif arises from dominant electrostatic O²⁻⋯X⁺ interactions embedded within a strongly covalent X–E scaffold that compensates intrinsic X–X repulsion, with only secondary multicenter O–X delocalization contributing to the overall electronic stabilization of the D_4h framework.

3.3 Magnetic response and assessment of aromaticity

To assess the magnetic response of the nine D_4h-ptO GMs, MICD analysis was performed (Fig. 5). The vector plots show that the induced currents do not exhibit a continuous global ring current; instead, two localized diatropic vortices are observed (Fig. 5a). Circuit i is confined to the region surrounding the central oxygen atom, whereas circuit ii is localized over the peripheral 3c X–E–X motifs of the crown (Fig. 5b). The corresponding ring current strengths (Fig. 5c) indicate that the peripheral circuit is dominant (10.5–13.0 nA T⁻¹ across the 3 × 3 series), while the oxygen-centered vortex is weaker (5.1–7.5 nA T⁻¹). This partitioning mirrors the bonding analysis: the covalent X–E framework (large WBI_(X–E), significant V^int_XC (X–E), and higher δ(X–E)) sustains the stronger local diatropic response, whereas the predominantly electrostatic O–X interactions (low WBI_(O–X), small δ (O–X), Coulomb-dominated IQA terms) support only a modest localized vortex. The MICD results, therefore, indicate superposition of local magnetic contributions rather than global aromatic stabilization of the D_4h framework.

Fig. 5 (a) MICD-vector plots of the O©X₄E₄²⁺ (X = Al–In, E = S–Te) systems. (b) Circuit detected of MICD-vector plots. (c) RCSs (in nA T⁻¹) for each circuit detected at the ωB97X-D/def2-TZVP//PBE0-GD3/def2-TZVP level.

Complementary to the MICD analysis, the global EDDB_G results (Fig. 6) further exclude global aromatic stabilization of the D_4h framework. For the representative sulfide series O©X₄S₄²⁺ (X = Al, Ga, In), the total populations of delocalized electrons are 5.983, 6.330, and 6.329 |e|, respectively, with the corresponding Se and Te analogues (Fig. S16) displaying comparable values and spatial distributions. The EDDB_G isosurfaces show that delocalized density is predominantly localized over the peripheral X–E scaffold, particularly along the 3c X–E–X connections, whereas only minimal contribution is found around the central oxygen. This distribution is fully consistent with the MICD results, which reveal dominant local diatropic vortices in the periphery and only a weaker localized response at the ptO center, and it corroborates the bonding analyses indicating that electron sharing is concentrated in the covalent X–E framework rather than in a globally delocalized O-centered circuit.

Fig. 6 Isosurface contour by global EDDB_G function (isovalue = ±0.006) with the corresponding population of delocalized electrons for O©Al₄S₄²⁺, O©Ga₄S₄²⁺ and O©In₄S₄²⁺ at the ωB97X-D/def2-TZVP level.

3.4 Exploring π-anion complexation as a stabilization strategy for ptO dications

An interesting extension of the present strategy would involve symmetric “full-sandwich” architectures in which the dicationic ptO core is framed by anions on both sides. However, the electrostatic potential (MEP) of the O©X₄E₄²⁺ clusters reveals a highly electropositive central region surrounded by comparatively more electronegative peripheral zones (Fig. 2b), so that the introduction of two anionic ligands leads to competing electrostatic effects, including anion–anion repulsion and disruption of the optimal electrostatic embedding that stabilizes the ptO motif. In contrast, the half–sandwich complexes examined here alleviate these repulsive interactions by allowing the external ligands to relax out of the plane and away from the π-surface of COT²⁻, thereby minimizing electrostatic conflict while preserving the planar X₄E₄ core. These observations suggest that viable full-sandwich architectures would require weakly coordinating counterions that maintain the balance between electrostatic confinement and structural rigidity.

Planar cyclooctatetraenide (COT²⁻), which adopts an aromatic D_8h-like geometry in its dianionic form, was selected as a representative π-delocalized counteranion. Upon complexation (Fig. 7), a marked dependence on the triel element is observed. For the aluminum derivatives, the central Al₄O core remains strictly planar: the O–Al distances (2.26 Å) and WBI(O–Al) values (0.20–0.21) are essentially identical to those of the isolated dications, and the charge distribution remains strongly polarized (q(O) ≈ −0.5; q(Al) = +1.2–1.5). Although the peripheral chalcogen atoms deviate from perfect coplanarity upon association, the square Al₄ framework and the ptO unit are preserved, indicating that the interaction with COT²⁻ is predominantly electrostatic and does not alter the intrinsic Al₄O bonding topology.

Fig. 7 Optimized structures of the O©X₄E₄²⁺-COT²⁻ adducts (X = Al, Ga, In; E = S–Te) and selected structural parameters (r(O–X), WBI(O–X), and atomic charges), illustrating preservation of the planar Al₄O core and distortion in the Ga and In analogues upon π-anion complexation.

In contrast, the gallium and indium derivatives exhibit substantial distortion of the square X₄ arrangement upon complexation. Despite moderate WBI(O–Ga) (0.49–0.52) and WBI(O–In) (0.42–0.51) values, the oxygen center no longer maintains an ideal planar tetracoordinate geometry. The comparatively lower positive charge at Ga and In (≈+1.0–+1.4) and the greater covalent contribution to O–X bonding correlate with increased structural flexibility and reduced resistance to external electrostatic perturbation.

Thus, π-anion complexation shows that retention of the ptO motif under charge compensation is not intrinsic to all systems but depends critically on the structural rigidity and electrostatic dominance of the Al₄O core. These findings indicate that counterion-assisted stabilization is feasible only when the cationic framework possesses sufficient internal robustness to preserve planar tetracoordination.

4 Conclusions

A systematic exploration of 25 O©X₄E₄²⁺ combinations (X = B, Al, Ga, In, Tl; E = O, S, Se, Te, Po) led to the identification of nine D_4h global minima featuring a genuine ptO center, specifically for X = Al, Ga, In and E = S, Se, Te. In these species, bonding is best described as a highly anionic O²⁻ unit electrostatically embedded within a charge-delocalized X₄E₄ tetracationic framework that is reinforced by a covalent X–E σ-scaffold, while direct X–X bonding remains negligible. Bonding analysis indicates that covalent electron sharing is concentrated within the peripheral three-center X–E–X units, whereas the oxygen center primarily retains localized lone pairs and engages in polarized O–X interactions. Magnetic response calculations reveal local diatropic vortices around both the oxygen center and the X–E–X regions; however, the absence of a continuous diatropic ring current across the D_4h framework excludes aromaticity.

Stabilization of planar tetracoordinate oxygen within the X₄E₄ platform is confined to a defined periodic regime. Only intermediate group-13 elements (Al–In) combined with mid-row chalcogens (S–Te) achieve the balance between electrostatic O²⁻ confinement and covalent σ-framework reinforcement required for a D_4h global minimum. Systems that are excessively covalent (B), insufficiently electrostatic (Tl), overly electronegative (O), or too diffuse (Po) fall outside this window. These results delineate clear compositional boundaries for electrostatic embedding and provide predictive criteria for future ptO design.

In addition to their thermodynamic stability as global minima, the dynamic stability of these systems is supported by BOMD simulations, which show structural persistence at finite temperatures. While kinetic pathways are not explicitly addressed, these results provide further support for their potential viability. Born–Oppenheimer molecular dynamics simulations demonstrate that all nine species remain structurally intact at 500 K, confirming that these ptO dications correspond to dynamically robust gas-phase minima. π-Anion complexation with COT²⁻ demonstrates that preservation of the planar ptO motif under external charge compensation is not intrinsic but depends critically on the rigidity and electrostatic dominance of the cationic framework, being maintained only for the aluminum derivatives. As a conceptual model for counterion confinement, these results indicate that ptO dications can, in principle, be trapped within charge-complementary environments provided that the internal scaffold is sufficiently robust to withstand electrostatic perturbation.

Collectively, these findings establish electrostatic embedding of O²⁻ within a covalently reinforced cationic X₄E₄ framework as a rational design strategy for stabilizing planar tetracoordinate oxygen without invoking classical covalent hypercoordination or global aromatic delocalization.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Figures S1–S15, Tables S1–S3 and cartesian coordinates. See DOI: https://doi.org/10.1039/d6ra04411e.

Acknowledgements

This work was supported by the financial support of the National Agency for Research and Development (ANID) through FONDECYT project 1241066 (W. T.) and Scholarship Program/BECAS DOCTORADO UNAB (D. S. S; L. M.; J. C.). This work was funded by the Postdoctoral Talent Attraction Competition for Research Centers and Institutes of the Universidad Andres Bello (UNAB) 2025, project N DI-02-25/ATP (L. L-P.). National Agency for Research and Development (ANID)/Scholarship Program/BECAS DOCTORADO NACIONAL/2024-21241825 (V. R-V). Powered@NLHPC: This research was partially supported by the NLHPC's supercomputing infrastructure (CCSS210001).

References

1. R. Hoffmann, R. W. Alder and C. F. Wilcox, J. Am. Chem. Soc., 1970, 92, 4992–4993.
2. L. Leyva-Parra, L. Diego, D. Inostroza, O. Yañez, R. Pumachagua-Huertas, J. Barroso, A. Vásquez-Espinal, G. Merino and W. Tiznado, Chem. Eur. J., 2021, 27, 16701–16706.
3. Z.-h. Cui, Y.-h. Ding, J. L. Cabellos, E. Osorio, R. Islas, A. Restrepo and G. Merino, Phys. Chem. Chem. Phys., 2015, 17, 8769–8775.
4. X. D. Jia and Z.-w. Du, Phys. Chem. Chem. Phys., 2023, 25, 4211–4215.
5. L. Leyva-Parra, D. Inostroza, O. Yañez, J. C. Cruz, J. Garza, V. García and W. Tiznado, Atoms, 2022, 10, 27.
6. M.-h. Wang, M. Orozco-Ic, L. Leyva-Parra, W. Tiznado, J. Barroso, Y. H. Ding, Z. H. Cui and G. Merino, J. Phys. Chem. A, 2021, 125, 3009–3014.
7. V. Vassilev-Galindo, S. Pan, K. J. Donald and G. Merino, Nat. Rev. Chem., 2018, 2, 0114.
8. R. Báez-Grez, A. Vásquez-Espinal and R. Pino-Rios, ChemPhysChem, 2025, 26, e202500350.
9. X. B. Liu, W. Tiznado, L. J. Cui, J. Barroso, L. Leyva-Parra, L. H. Miao, H. Y. Zhang, S. Pan, G. Merino and Z. H. Cui, J. Am. Chem. Soc., 2024, 146, 16689–16697.
10. L. Leyva-Parra, L. Diego, O. Yañez, D. Inostroza, J. Barroso, A. Vásquez-Espinal, G. Merino and W. Tiznado, Angew. Chem., Int. Ed., 2021, 60, 8700–8704.
11. J. K. Burdett and T. Hughbanks, J. Am. Chem. Soc., 1984, 106, 3101–3113.
12. S. D. Li, G. M. Ren, C. Q. Miao and Z. H. Jin, Angew. Chem., Int. Ed., 2004, 43, 1371–1373.
13. S. D. Li, C. Q. Miao and G. M. Ren, Eur. J. Inorg. Chem., 2004, 2232–2234.
14. K. Zhou, S. Roy and C. B. Zhao, Comput. Theor. Chem., 2017, 1110, 35–39.
15. K. Zhou, S. K. Roy and C. B. Zhao, Russ. J. Inorg. Chem, 2019, 64, 303–307.
16. G. Castillo-Toraya, F. Ortíz-Chi, J. Barroso, M. Orozco-Ic, L. Leyva-Parra and G. Merino, Angew. Chem., Int. Ed., 2025, 64, e202500292.
17. B. Jin, T. Li, M. Yan and Y. J. Wang, J. Phys. Chem. A, 2026, 130, 830–835.
18. R. Sun, Y. Yang, X. Wu, H. J. Zhai, C. Yuan and Y.-B. Wu, Chem. Sci., 2025, 16, 12873–12878.
19. K. Franzreb and P. Williams, Chem. Phys. Lett., 2006, 419, 379–384.
20. R. Sun, C. Yuan, H. J. Zhai and Y.-B. Wu, J. Chem. Phys., 2023, 158, 144301.
21. X. D. Jia and Z. W. Du, Phys. Chem. Chem. Phys., 2024, 26, 3907–3911.
22. D. S. Sacanamboy, V. Roman-Ventura, O. Yañez, W. García-Argote, L. Leyva-Parra and W. Tiznado, Phys. Chem. Chem. Phys., 2026, 28, 7277–7281.
23. O. Yañez, R. Báez-Grez, D. Inostroza, W. A. Rabanal-León, R. Pino-Rios, J. Garza and W. Tiznado, J. Chem. Theory Comput., 2019, 15, 1463–1475.
24. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
25. W. Küchle, M. Dolg, H. Stoll and H. Preuss, Mol. Phys., 1991, 74, 1245–1263.
26. W. Küchle, M. Dolg, H. Stoll and H. Preuss, J. Chem. Phys., 1994, 100, 7535–7542.
27. P. Fuentealba, H. Preuss, H. Stoll and L. Von Szentpály, Chem. Phys. Lett., 1982, 89, 418–422.
28. P. Fuentealba, L. Von Szentpaly, H. Preuss and H. Stoll, J. Phys. B: At., Mol. Opt. Phys., 1985, 18, 1287.
29. A. Bergner, M. Dolg, W. Küchle, H. Stoll and H. Preuß, Mol. Phys., 1993, 80, 1431–1441.
30. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132.
31. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
32. C. Riplinger and F. Neese, J. Chem. Phys., 2013, 138, 034106.
33. C. Riplinger, B. Sandhoefer, A. Hansen and F. Neese, J. Chem. Phys., 2013, 139, 134101.
34. C. Riplinger, P. Pinski, U. Becker, E. F. Valeev and F. Neese, J. Chem. Phys., 2016, 144, 024109.
35. D. G. Truhlar, Chem. Phys. Lett., 1998, 294, 45–48.
36. F. Neese, A. Hansen and D. G. Liakos, J. Chem. Phys., 2009, 131, 064103.
37. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson and H. Nakatsuji, Gaussian 16, Revision C.02, Gaussian, Inc., Wallingford, CT, 2016.
38. F. Neese, WIREs Comput Mol Sci., 2022, 12, e1606.
39. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., 1989, 36, 199–207.
40. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735–746.
41. K. B. Wiberg, Tetrahedron, 1968, 24, 1083–1096.
42. NBO 7.0. E. D. Glendening, J. , K. Badenhoop, A. E. Reed, J. E. Carpenter, J. A. Bohmann, C. M. Morales, P. Karafiloglou, C. R. Landis, and F. Weinhold, Theoretical Chemistry Institute, University of Wisconsin, Madison ( 2018).
43. E. D. Glendening, C. R. Landis and F. Weinhold, J. Comput. Chem., 2019, 40, 2234–2241.
44. D. Y. Zubarev and A. I. Boldyrev, Phys. Chem. Chem. Phys., 2008, 10, 5207–5217.
45. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
46. A. M. Pendas, M. A. Blanco and E. Francisco, J. Chem. Phys., 2004, 120, 4581–4592.
47. M. A. Blanco, A. Martín Pendás and E. Francisco, J. Chem. Theory Comput., 2005, 1, 1096–1109.
48. A. M. Pendás, M. A. Blanco and E. Francisco, J. Comput. Chem., 2007, 28, 161–184.
49. AIMAll (Version 19.10.12), T. A. Keith, TK Gristmill Software, Overland Park KS, USA, 2019.
50. T. A. Keith, Chem. Phys., 1996, 213, 123–132.
51. T. A. Keith and R. F. W. Bader, Chem. Phys. Lett., 1993, 210, 223–231.
52. T. A. Keith and R. F. W. Bader, Chem. Phys. Lett., 1992, 194, 1–8.
53. D. Sundholm, H. Fliegl and R. J. F. Berger, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2016, 6, 639–678.
54. H. Fliegl, S. Taubert, O. Lehtonen and D. Sundholm, Phys. Chem. Chem. Phys., 2011, 13, 20500–20518.
55. J. Jusélius, D. Sundholm and J. Gauss, J. Chem. Phys., 2004, 121, 3952–3963.
56. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
57. D. W. Szczepanik, M. Solà, M. Andrzejak, B. Pawełek, J. Dominikowska, M. Kukułka, K. Dyduch, T. M. Krygowski and H. Szatylowicz, J. Comput. Chem., 2017, 38, 1640–1654.
58. K. Wolinski, J. F. Hinton and P. Pulay, J. Am. Chem. Soc., 1990, 112, 8251–8260.
59. D. W. Szczepanik, M. Andrzejak, K. Dyduch, E. Żak, M. Makowski, G. Mazur and J. Mrozek, Phys. Chem. Chem. Phys., 2014, 16, 20514–20523.
60. D. W. Szczepanik, M. Andrzejak, J. Dominikowska, B. Pawełek, T. M. Krygowski, H. Szatylowicz and M. Solà, Phys. Chem. Chem. Phys., 2017, 19, 28970–28981.

---