A ternary Li₃B₇Si₃ cluster: when reduced electron delocalization enhances thermodynamic stability

Abstract

The stability and electronic structure of B₇Si₃^q (q = +, 0, −, 2− and 3−), Li₃B₇Si₃ and Li₃B₁₀H₃ clusters have been investigated through a comprehensive isomer search, with a focus on the perfect triangular planar isomer T⁺ of B₇Si₃⁺. Despite its double aromaticity, this form is thermodynamically unstable. Two strategies for the stabilization of the triangular form are proposed: (i) sequential addition of σ electrons and (ii) incorporation of Li atoms to reduce excessive negative charge. The former reveals that the addition of up to four electrons gradually improves its stability, with T²⁻ being the most favourable isomer in this charged state. The latter strategy leads to the star form emerging as a global minimum of the ternary Li₃B₇Si₃ cluster, which is significantly more stable. To elucidate the origin of these stabilizations, ring current maps and AdNDP, ELF_σ and ETS-NOCV analyses were employed. Results show a shift from the dominant σ aromaticity in T⁺ to a more balanced σ–π aromatic character in Li₃B₇Si₃, akin to benzene. The AdNDP analysis reveals that T⁺ possesses four σ delocalized bonds, whereas Li₃B₇Si₃ contains only three, serving as one of the indicators of reduced σ-delocalization. ELF_σ analysis further revealed an enhanced peripheral σ delocalization in Li₃B₇Si₃, highlighting the critical role of edge-localized electrons in stabilizing planar aromatic clusters.

---

1. Introduction

As one of the simplest aromatic hydrocarbons, benzene holds a unique position in the field of chemistry. Since its discovery in 1825 by Michael Faraday,^1,2 benzene has continuously challenged scientists due to its unusual geometry, electronic structure, high thermodynamic stability, and peculiar chemical behavior.^3–5 It took more than a century from that initial discovery of benzene to the formulation by Hückel of the electron counting rule of (4N + 2) in the early 1930s,^6–8 which has since provided chemists with a theoretical framework to rationalize its high stability through the concept of aromaticity defined by cyclic delocalization of π electrons. This theoretical milestone has fundamentally advanced the understanding of the chemical bonding phenomena and aromatic behavior of organic compounds. Nowadays, the aromatic character has also been extensively applied to different classes of atomic clusters that have been assigned not only as a single π aromatic, but also as a single σ aromatic,^9–14 double (σ + π) aromatic, or even as a δ aromatic.^15–17

Several factors contribute to the stabilization of a molecular system. It is generally observed that in atomic clusters where aromaticity is present, such as in benzene, there is a dominant stabilizing factor; most of the reported stable clusters exhibit a double aromaticity, and this characteristic is often regarded as a hallmark of their structural stability. In contrast, less stable structures are typically characterized by a lack of electron delocalization or, in other words, by an absence of aromaticity. In some rare cases, a cluster isomer that exhibits double (σ + π) aromaticity is found to be unstable. This raises a legitimate question as to whether a decrease in electron delocalization, of either the σ or π system, in a structure is associated with enhanced thermodynamic stability. This is an apparently absurd question within the chemical orthodoxy, but it is worth raising. In this context, we set out to search for such a compound using quantum chemical computations. Our extensive search led to the stable configurations of the binary B₇Si₃⁺ cluster cation and subsequently, the ternary Li₃B₇Si₃ cluster. Herein, we report a strategy for stabilizing this type of isomer and explore how variations in the σ electron delocalization influence the stability of planar clusters.

2. Computational methods

Each cluster considered is individually explored through a combination of a stochastic algorithm¹⁸ and a genetic algorithm¹⁹ for structural generation. Initial geometries are optimized using density functional theory with the PBE²⁰ functional and the 6-31G(d) basis set,^21–23 without performing vibrational frequency calculations. Stationary points located within 3 eV from the most stable structure are further re-optimized using the same PBE functional but with the larger 6-311+G(d) basis set,^21–23 and harmonic vibrational analyses are subsequently carried out to identify the energy minima structures. The most stable isomers are validated by single-point electronic energy, with coupled-cluster theory CCSD(T) computations using the aug-cc-pVTZ basis set. The Gaussian 16 program²⁴ is used to carry out all relevant electronic structure calculations.

Using the BHandHLYP functional with the 6-311+G(d) basis set, origin-independent current densities CTOCD-DZ2 (ref. 25) and bond current strengths (BCS)²⁶ were calculated, and these parameters were further decomposed into π and σ contributions. The BHandHLYP functional was calibrated by a recent benchmark study for reliably describing magnetic properties.²⁷ Magnetic current density maps and BCS were computed and visualized by using the SYSMOIC²⁸ software package. All current density maps of planar structures were computed in the molecular plane and in a plane located at 1 Bohr above and parallel to the molecular plane, under the influence of a magnetic field applied perpendicular to the molecular plane. The ring current is characterized by following the convention on the direction of the induced current, that is, a clockwise circulation indicates a diatropic current, which is associated with aromatic character, whereas a counterclockwise circulation corresponds to a paratropic current, which is associated with antiaromatic character.

Chemical bonding characteristics are also analyzed using the Adaptive Natural Density Partitioning (AdNDP) method²⁹ as implemented in the Multiwfn program suite.^30,31 Electron localization function (ELF)^32,33 maps are constructed using the Dgrid 5.0 package³⁴ and visualized by ChimeraX software.³⁵

Energy decomposition analysis (EDA) combined with the natural orbital for the chemical valence (NOCV) method was performed at the B3LYP-D3BJ/def2-QZVP³⁶ level using the ORCA program package.³⁷ In this scheme, the interaction energy (ΔE_int) is partitioned into electrostatic (ΔE_elstat), Pauli repulsion (ΔE_Pauli), orbital (ΔE_orb), and dispersion (ΔE_disp) contributions.

Ab initio molecular dynamics (AIMD) simulations are carried out with the ORCA program package,³⁷ using the PBE exchange–correlation functional with D3 dispersion correction and the def2-TZVPD basis set.^38,39 The simulations were performed using a Nosé–Hoover chain thermostat at 300 K. A time step of 2.0 fs was employed to simulate 20 ps of AIMD to 20 ps, with atomic positions recorded every 10 fs.

3. Results and discussion

3.1. Instability of a doubly aromatic cluster structure

Geometrical structures of the most stable isomers of the binary clusters B₇Si₃^q in five different charge states, q = +, 0, −, 2−, and 3−, are presented in Fig. 1 with additional structures provided in Fig. S1. The isomers are labeled B₇Si₃^q.X with X = 1, 2, 3, … to denote their ordering in relative energy with respect to the most stable structure. The global energy minimum B₇Si₃⁺.1 of the cation is formed by the attachment of two Si atoms to the edge of a B₇Si umbrella-like structure.⁴⁰ Accordingly, B₇Si₃⁺.1 is more stable than the second lowest-lying isomer B₇Si₃⁺.2 by up to 14 kcal mol⁻¹ in relative energy, making it particularly noteworthy for an in-depth analysis. However, we focus in the present study on a different type of isomer denoted as T^q (T denotes a triangle shape) in Fig. 1. Let us consider the cation T⁺, which exhibits a perfect triangular planar shape resulting from the attachment of three Si atoms to three opposite edges of a B₇ unit. It possesses an imaginary frequency of 243i cm⁻¹ and is less stable by ∼40 kcal mol⁻¹ than the corresponding lowest-lying B₇Si₃⁺.1 isomer.

Fig. 1 Optimized geometries of the B₇Si₃^q cluster (q = +, 0, −, 2−, and 3−) obtained at the PBE/6-311+G(d) level. Relative energies (kcal mol⁻¹) were calculated using (U)CCSD(T)/aug-cc-pVTZ//PBE/6-311+G(d), with values from PBE/6-311+G(d) given in parentheses. All energies include zero-point energy (ZPE) corrections evaluated at the PBE/6-311+G(d) level.

Fig. 2 presents three π-MOs and frontier σ-MOs, including the doubly degenerate σ LUMO of T⁺ (B₇Si₃⁺.13). The addition of electrons to T⁺ to transform the degenerate LUMO into the HOMO can increase only the number of σ electrons but retains its 6 π electrons. Consequently, this allows us to figure out how changes in the σ electron number can influence structural stability. To further investigate this issue, we incrementally added electrons to the B₇Si₃⁺. Stable structures of B₇Si₃^q are also presented in Fig. 1. The addition of one electron to T⁺ eliminates the imaginary frequency, making the corresponding neutral T (B₇Si₃.11) a local energy minimum and reducing its relative energy with respect to the global minimum B₇Si₃.1 to 19 kcal mol⁻¹. Further electron additions continue to reduce the relative energy to ∼12 and ∼4 kcal mol⁻¹ for T⁻ (B₇Si₃⁻.7) and T²⁻ (B₇Si₃²⁻.3), respectively, regarding their global minima. Upon full occupation of the degenerate LUMO of T⁺ with four electrons, the resulting trianion T³⁻ (B₇Si₃³⁻.4) becomes less stable than B₇Si₃³⁻.1 by ∼4 kcal mol⁻¹, along with the emergence of double-chain ribbon isomers such as B₇Si₃³⁻.1 and B₇Si₃³⁻.3. This observation is in line with earlier findings suggesting that the accumulation of electrons in an atomic cluster generally leads to a reduction in structural dimensionality; that is, going from 3D → 2D → 1D.⁴¹

Fig. 2 Three π-MOs and the frontier σ-MO shapes of T⁺.

In view of the inherent instability of highly charged polyanions with respect to electron rejection, a commonly employed approach to stabilize them is the incorporation of alkali metals into the structures. As shown in Fig. 3, this approach proves to be effective, with three lithium atoms supplying three electrons to the cluster.

Fig. 3 Optimized geometries of Li₃B₇Si₃ (abbreviated as LiSi) and Li₃B₁₀H₃ (abbreviated as LiH) obtained at the PBE/6-311+G(d) level. Relative energies (kcal mol⁻¹) were calculated using (U)CCSD(T)/aug-cc-pVTZ//PBE/6-311+G(d), with values from PBE/6-311+G(d) given in parentheses. All energies include zero-point energy (ZPE) corrections evaluated at the PBE/6-311+G(d) level. The natural atomic charges of LiSi.1 and LiH.1 are also provided.

An isolobal substitution Si = [BH] is thus carried out, revealing perfectly planar D_3h LiSi.1 (Li₃B₇Si₃ is abbreviated as LiSi) and LiH.1 (Li₃B₁₀H₃ is abbreviated as LiH) structures emerging as the global energy minima for both Li₃B₇Si₃ and Li₃B₁₀H₃, respectively. Such planar arrangements are exceptionally rare because Li⁺ ions typically favour the formation of three-dimensional structures when bonding to boron clusters.^42,43 Moreover, a planar geometry of boron clusters is advantageous for their synthesis via deposition onto solid surfaces.^44,45 The substitution of three Si atoms with three [BH] units increases the number of atoms in the cluster by three, resulting in a greater number of possible isomers. Nevertheless, Li₃B₁₀H₃ still prefers a planar configuration as the most stable form, LiH.1, with the next-lowest-energy isomer LiH.2 adopting a distorted geometry and lying only ∼4 kcal mol⁻¹ above LiH.1. In this study, we focus on elucidating the stability of LiSi.1, while the case of LiH.1 is explained in an entirely analogous manner, with the corresponding discussion provided in the SI.

The T₁ diagnostic values from single-point CCSD(T)/aug-cc-pVTZ calculations for the structures presented in Fig. S1 and 3 are summarized in Table 1. The results indicate that several isomers exhibit T₁ values exceeding the thresholds of 0.02 for closed-shell and 0.04 for open-shell systems.⁴⁶ However, the most stable isomers (B₇Si₃⁺.1, B₇Si₃.1, B₇Si₃⁻.1, B₇Si₃²⁻.1, and B₇Si₃³⁻.1), as well as the triangular form (T^q), display T₁ values within the reliability limits for a single-reference description. Notably, upon the addition of three Li atoms, the structures exhibited an even more pronounced single-reference character, as reflected by the lower T₁ values of Li₃B₇Si₃ and Li₃B₁₀H₃ isomers.

Table 1 Bond current strengths (BCS) for selected bonds. The values are presented separately for σ and π electron contributions, including their relative ratios, as well as for the total contributions from all electrons. All BCS values are given in nA T⁻¹

Species	Electron	B1–B2	B2–B3	B1–Si1	Total of B ring^a
a The B ring is the B1B2B3B4B5B6 ring.
T^+	σ	14.9	17.5	7.0	97.2
	π	3.7	4.8	0.6	25.5
	σ/π	4.0	3.6	11.7	3.8
	All	18.6	22.3	7.6	122.7
LiSi.1	σ	5.2	14.9	5.2	60.3
	π	4.7	5.8	1.2	31.5
	σ/π	1.1	2.6	4.3	1.9
	All	9.9	20.7	6.4	91.8
LiH.1	σ	5.9	14.0	10.7	59.7
	π	5.3	6.4	2.3	35.1
	σ/π	1.1	2.2	4.7	1.7
	All	11.2	20.4	12.7	94.8

---

3.2. Chemical bonding and aromaticity

The most fundamental difference between T⁺ and LiSi.1 is that the doubly degenerate σ LUMO of T⁺ becomes the doubly degenerate σ HOMO of LiSi.1, as shown in Fig. 4. The following analysis aims to clarify the effect of the addition of σ electrons to the structure, which transforms the LUMO into the HOMO.

Fig. 4 The HOMO, LUMO and three π-MOs shapes of LiSi.1.

Firstly, the magnetic current density maps for T⁺ and LiSi.1 are presented in Fig. 5, in comparison to those of benzene, showing contributions from σ, π, and all electrons at 1 Bohr above the molecular plane, and the maps contributed from σ electrons located in molecular planes are shown in Fig. S2 (SI). Their bond current strengths (BCS) are shown in Fig. 6; each arrow in the figure represents the net current strength as a percentage relative to the reference value of 12 nA T⁻¹, which corresponds to the net diatropic current strength of benzene.⁴⁷ The corresponding figures for Li₃B₁₀H₃.1 are also provided in the SI as Fig. S3 and S4. The BCS of selected bonds for T⁺, LiSi.1, and Li₃B₁₀H₃.1 are also presented in Table 1, with units in nA T⁻¹.

Fig. 5 Magnetic current density maps of T⁺ (top), LiSi.1 (middle) and benzene (bottom) from σ-electron contribution (left), π-electron contribution (center), and all-electron contribution (right) at 1 Bohr above the molecular planes. Diatropic and paratropic domains are marked by green and red rings, respectively.

Fig. 6 Bond current strengths of T⁺ (top), LiSi.1 (middle) and benzene (bottom) from σ electron contribution (left), π electron contribution (center), and all electrons (right). Each arrow represents the percentage relationship with respect to the diatropic benzene current strength of 12 nA T⁻¹.

From the magnetic ring current maps of T⁺ shown in Fig. 5, the first remarkable observation is that the T⁺ cation exhibits double aromaticity, even though it possesses an imaginary frequency and a high energy location. The aromatic character of different compounds can be compared by evaluating the magnitude of BCS values calculated for the bonds in the B1B2B3B4B5B6 ring of T⁺ and LiSi.1 (cf. Fig. 6 for atom labels), relative to the C–C bonds in benzene. Ring current maps indicate that the σ aromaticity is significantly pronounced, whereas the π-aromaticity is relatively weak. Although T⁺ possesses, like benzene, 6 π-electrons, the π BCSs for the B1–B2 and B2–B3 bonds of T⁺ amount to 3.7 and 4.8 nA T⁻¹, respectively, whereas π electrons in benzene contribute 5.8 nA T⁻¹ to C–C bonds. While the σ electrons in benzene are localized, the diatropic domain in the outer ring is larger than the paratropic domain in the inner ring, resulting in large σ BCSs for C–C bonds, reaching a value of 6.1 nA T⁻¹.²⁸ The σ electron ring current map of T⁺ shows a strong diatropic current throughout the entire structure, without any paratropic domain observed, indicating that σ electrons in this structure are mostly delocalized. The BCS values of the B1–B2 and B2–B3 bonds formed from σ electrons in T⁺ are up to remarkably high values of 14.9 (126%) and 17.5 nA T⁻¹ (147%), respectively.

The total BCS of the whole B1B2B3B4B5B6 ring in T⁺ of σ electrons of 123 nA T⁻¹ is nearly four times as great as that of the corresponding π electrons (97 nA T⁻¹). Thus, an imbalance emerges between delocalization of the σ and π electrons in the T⁺ species. If this imbalance can be reduced, leading to an increased structural stability, it can be viewed that a planar geometry can confer a greater thermodynamic stability when the aromatic contributions from σ and π electrons become more balanced.

A similar analysis was conducted for LiSi.1, in which four σ electrons were added from the electron transfer of Li atoms to the two-fold degenerate LUMO (5e′) of T⁺ (cf. Fig. 2), turning it into the two-fold degenerate HOMO in LiSi.1 (cf. Fig. 4), revealing a shift in behaviour that more closely resembles that of benzene. As six π electrons remain, the π aromaticity in LiSi.1 slightly increases compared to that in T⁺, as reflected by the π-BCS values for the B1–B2 and B2–B3 bonds that amount to 4.7 and 5.8 nA T⁻¹, respectively (cf. Table 1). In contrast, despite an increased number of σ electrons, the σ-BCS values for the B1–B2 and B2–B3 bonds tend to decrease to 5.2 and 14.9 nA T⁻¹, respectively. The reduction in σ electron BCS values at the B1–B2 bond is substantial in going from 14.9 nA T⁻¹ in T⁺ to 5.2 nA T⁻¹ in LiSi.1.

While a strong diatropic outer ring current near B1–B2 can still be observed in LiSi.1, the emergence of a significant paratropic inner ring current (indicated by the red circle in Fig. 5) contributes largely to this BCS decrease. The B atom in the [BH] group appears to facilitate a better π-electron delocalization from its 2p_z orbital in LiH.1 than contributions from the 3p_z orbital of Si and the 2p_z orbital of B in LiSi.1. As a result, the total BCS value of LiH.1 (35.1 nA T⁻¹) is higher than that of LiSi.1 (31.5 nA T⁻¹). The BCS contribution from σ electrons in LiH.1 (59.7 nA T⁻¹) was observed to be slightly lower than that in LiSi.1 (60.3). Overall, these results indicate that the σ BCS to π BCS ratio (cf. Table 1) in LiH.1 (∼1.7) is smaller than that in LiSi.1 (∼1.9). These preliminary results suggest that LiH.1 possesses a higher thermodynamic stability than LiSi.1, which will be further verified by the ELF analysis presented in the following section.

This observation also suggests a shift in the σ electron delocalization from the entire structure in T⁺ toward the periphery in LiSi.1. This hypothesis is further supported by AdNDP and ELF_σ analyses presented in the following paragraphs.

The AdNDP analyses of T⁺ and LiSi.1 are presented in Fig. 7, and LiH.1 is presented in Fig. S5. Each structure possesses six π electrons, thus satisfying Hückel's counting rule. Their aromatic character is clearly manifested, as illustrated in the corresponding ring current maps (cf. Fig. 5). As initially expected, changes in the number of σ electrons significantly affect the structural stability. The lone pairs on Si and localized electrons in both T⁺ (cf. (1), (2), and (3) bonds) and LiSi.1 (cf. (7), (8) and (9) bonds) are quite similar to each other, with bonds in LiSi.1 containing higher occupation numbers (ONs). T⁺ contains one 9c–2e σ bond (5) and three 10c–2e σ bonds (6), indicating that T⁺ possesses eight delocalized σ electrons. Upon the addition of four σ electrons to T⁺, the resulting LiSi.1 structure did not contain 12 delocalized σ electrons as might be expected, but instead only six delocalized (11c–2e) σ electrons as indicated in bonds (11). The remaining six localized (3c–2e) σ electrons of LiSi.1 formed three σ bonds as indicated in bonds (10) (cf. Fig. 7). This suggests that the addition of σ electrons to T⁺ induces a radical reorganization of its σ electronic structure, redistributing the delocalized and localized electrons. As a result, LiSi.1 exhibits fewer delocalized σ and more localized σ electrons, which is reflected in the change in aromatic character observed in the B1B2B3B4B5B6 ring. The results of the AdNDP analysis of LiH.1 yield are internally consistent with those of LiSi.1, demonstrating the isolobal relationship between these two isomers.

Fig. 7 AdNDP bonding patterns for (a) the T⁺ isomer of B₇Si₃⁺ and (b) the global minimum LiSi.1 of Li₃B₇Si₃. The occupation numbers (ONs) are indicated.

The change in the behaviour of σ electrons in going from T⁺ to LiSi.1 is remarkable. Therefore, analysis of the electron localization function of σ electrons (ELF_σ) was carried out for these two species to further elucidate this specific behaviour. The ELF_σ for both T⁺, LiSi.1, and LiH.1 are presented in Fig. 8 at different isosurface values, with red circles (with the label from bif(a) to bif(h)) added to highlight the locations of bifurcations corresponding to each isosurface value.

Fig. 8 Electron localization function ELF_σ maps of T⁺ (top), LiSi.1 (middle), and LiH.1 (bottom) at different isosurface values. Red circles (with the label from bif(a) to bif(h)) indicate the positions of bifurcations corresponding to each isosurface value.

The determining bifurcation values for both T⁺ and LiSi.1 are all greater than 0.8, indicating a high degree of delocalization of σ electrons within these structures.⁴⁸ The bifurcation observed at the bif(a) positions reflects the electronic separation of either the Si lone pair (in T⁺ and LiSi.1) or the B–H 2c–2e bond (in LiH.1) from the rest of the molecular domain. The bifurcation observed at the bif(c) positions reflects delocalization in peripheral regions. The bifurcations at other positions correspond to a delocalization inside the structure. Except for the bif(a) at the critical ELF = 0.63 in LiH.1, the other domains split at critical ELF values above 0.8, indicating a high degree of electron delocalization. The highest bifurcation value of T⁺ is 0.91 at the bif(d) position. The bif(d) position corresponds to the high degree of electron delocalization observed in the delocalized bond (5), as revealed by the AdNDP analysis of T⁺ (Fig. 7).

A similar delocalization bond (5) is not present in LiSi.1; instead, it is replaced by three localized 3c–2e bonds (bond (10) in Fig. 7), indicating a reduction to a lower critical ELF value of 0.85 at bif(f) of LiSi.1. This transformation also leads to the formation of three localized 3c–2e basins associated with the B1B2B7, B3B4B7, and B5B6B7 triangles in the LiSi.1 and LiH.1 structures, resulting in a more uniform distribution of σ-electrons throughout the framework. This observation suggests that excessive electron delocalization can reduce structural stability, and the redistribution of electrons to lower aromatic character can be an effective means of enhancing stability. Thus, although the degree of σ aromaticity tends to reduce in going from T⁺ to LiSi.1 (and LiH.1), as reflected by a marginal reduction in the highest bifurcation value of ELF_σ from 0.91 to 0.89 (and 0.88), the structural stability increases.

To gain deeper insights into the bonding characteristics, ETS-NOCV analysis was carried out for LiSi.1 using the [Li₃]³⁺ and [B₇Si₃]³⁻ fragmentation scheme, based on the natural atomic charges from the NBO analysis, as illustrated in Fig. 3. The results, as summarized in Table 2, reveal that the interaction is strongly dominated by the electrostatic attraction (ΔE_elstat = −704 kcal mol⁻¹), consistent with the large opposite charges of the fragments. In addition, a substantial orbital contribution (ΔE_orb = −185 kcal mol⁻¹) indicates significant orbital overlap and charge-transfer effects, highlighting that the bonding is not purely ionic but involves a noticeable degree of covalency. The Pauli repulsion is relatively small (+76 kcal mol⁻¹) compared to the attractive terms, and dispersion plays only a minor stabilizing role (−10 kcal mol⁻¹). Overall, the total interaction energy (ΔE_int = −823 kcal mol⁻¹) points to an exceptionally strong and stable bonding situation that combines predominant ionic character with appreciable covalent stabilization. The first five deformation densities and the corresponding eigenvalues and energetic contributions for LiSi.1 are depicted in Fig. S6 (SI). Among them, the first three σ-type channels exhibit very large eigenvalues (∼0.92 e) with a combined stabilization energy of −308 kcal mol⁻¹, indicating that σ interactions dominate the bonding. The deformation density plots (Δρ₃) reveal charge depletion (red) near the Li centers and charge accumulation (blue) in the vicinity of Si, consistent with a pronounced σ back-donation from the Li atoms directly into the orbitals of Si. In contrast, the two π-type channels have smaller eigenvalues (∼0.24 e) and minor energetic contributions (∼6–7 kcal mol⁻¹ each), suggesting only a secondary role in the overall interaction. These results clearly demonstrate that the [Li₃–B₇Si₃] bonding is primarily governed by strong electrostatic attraction, complemented by significant σ back-donation, while π contributions remain negligible. Importantly, this finding is consistent with the above conclusion regarding the attenuation of σ aromaticity, which in turn enhances the overall stability of LiSi.1.

Table 2 Results of EDA (B3LYP-D3(BJ)/def2-QZVP) for LiSi.1 with [Li₃]³⁺ + [B₇Si₃]³⁻ as fragments. Energy values are in kcal mol⁻¹

ΔEint	ΔEorb	ΔEelstat	ΔEPauli	ΔEdisp
−822.7	−184.9	−704.4	76.1	−10.4

---

The potential energy profile of LiSi.1 at 300 K (Fig. 9) exhibits stable fluctuations around a constant mean value without significant drifts or abrupt changes, confirming that the planar structure remains dynamically robust throughout the 20 ps AIMD trajectory. In addition, AIMD trajectory analyses from both top and side views (see SI: file Li3B7Si3.mp4) clearly reveal that the Li atoms exhibit small oscillatory excursions above and below the boron plane, but no net out-of-plane migration was observed; the oscillations are reversible and the Li atoms returned to their equilibrium positions throughout the 20 ps simulation.

Fig. 9 Potential energy trajectories of LiSi.1 at 300 K.

4. Concluding remarks

In the present theoretical study, we first conducted a structural isomer search to identify the stable configurations of the B₇Si₃^{+/0/−/2−/3−} in five charge states, Li₃B₇Si₃ and Li₃B₁₀H₃ clusters. While the survey revealed that the most stable isomer of B₇Si₃⁺ has a ribbon-type shape, our study focuses on the T⁺ isomer of B₇Si₃⁺, whose star shape exhibits properties contrary to conventional expectations. Specifically, this T⁺ cationic isomer contains a double aromaticity but turns out to be an unstable structure.

Given that the LUMO of the T⁺ isomer is a doubly degenerate σ MO, the addition of four electrons is considered an effective approach to stabilize this star form. Incremental electron addition was carried out up to the trianionic charge state, and the relative energy between 1^q and the lowest-lying isomer of B₇Si₃^q was found to consistently decrease in going to the dianion T²⁻ and the trianion T³⁻.

A second strategy for stabilization involves the replacement of the polyanionic charges by introducing lithium atoms that are capable of donating electrons while reducing the overall negative charge. This approach leads to the ternary Li₃B₇Si₃ cluster whose global energy minimum remains planar, given that the T⁺ form is originally less stable than the lowest-energy isomer B₇Si₃⁺.1 by as much as 40 kcal mol⁻¹.

In employing five analytical methods, namely the ring current map, bond current strength (BCS), AdNDP, ELF, and ETS-NOCV, the origin of the enhanced stability resulting from the addition of σ electrons was determined. The ring current maps point out that the LiSi.1 species retains a double (σ + π) aromaticity, similar to that of T⁺. However, a key difference is observed: in T⁺, the σ aromaticity is predominant over the π aromaticity, while in LiSi.1, the π aromaticity is enhanced while the σ aromaticity is reduced, approaching a BCS balance between both sets of σ and π electrons, as seen in benzene.

The AdNDP analysis indicates that the T⁺ form possesses four delocalized σ bonds, a feature that contributes to its instability. Although the electronic structure of LiSi.1 arises from the addition of four σ electrons to the doubly degenerate LUMO of T⁺, the AdNDP results show that LiSi.1 features only three delocalized σ MOs, consistent with the stability expected from Hückel's counting rule. This suggests that the resulting structure not only gains electrons but also undergoes a reorganization of its σ electron framework. Subsequent ELF_σ analysis is consistent with all other analytical methods, indicating that a reduction in σ electron delocalization contributes to an enhanced structural stability.

The 20 ps AIMD simulations at 300 K demonstrate that LiSi.1 retains a dynamically robust planar structure, with Li atoms exhibiting only small reversible oscillations above and below the boron plane without any net out-of-plane migration.

Conflicts of interest

There is no conflict of interest to declare.

Data availability

The quantum chemical program Gaussian 16 has been used. Reference of this program is given in the list of references.

The additional datasets supporting this article are given as part of the SI: Supplementary information: Cartesian coordinates of the optimized geometries, information on the electronic structures, magnetic ring current maps, electron localization function maps, and NOCV deformation densities of the clusters considered. See DOI: https://doi.org/10.1039/d5ra04954g.

Acknowledgements

MPPH acknowledges Ho Chi Minh City University of Technology HCMUT, Vietnam National University VNU-HCM, for supporting this study. MTN is indebted to Vin University for a Distinguished Professorship under grant no. 437/2025. The computer resources were supported by the PL-GRID infrastructure.

References

1. E. Mitscherlich, Ann. Pharm., 1834, 9, 39–48.
2. M. Faraday, Philos. Trans. R. Soc. London, 1825, 115, 440–466.
3. A. Kekulé, Bulletin de la Société chimique de Paris, 1865, 3, 98–110.
4. A. Ladenburg, Berichte der Deutschen Chemischen Gesellschaft, 1869, 2, 140–142.
5. J. Dewar, Proc. R. Soc. Edinburgh, 1869, 6, 82–86.
6. E. Hückel, Z. Med. Phys., 1931, 70, 204–286.
7. G. Frenking, Theor. Chem. Acc., 2000, 187–189.
8. J. A. Berson, Chemical Creativity: Ideas from the Work of Woodward, Hückel, Meerwein, and Others, Wiley, 1999.
9. E. Jimenez-Izal and A. N. Alexandrova, Phys. Chem. Chem. Phys., 2016, 18, 11644–11652.
10. R. W. A. Havenith, F. De Proft, P. W. Fowler and P. Geerlings, Chem. Phys. Lett., 2005, 407, 391–396.
11. T. J. Robilotto, J. Bacsa, T. G. Gray and J. P. Sadighi, Angew. Chem., Int. Ed., 2012, 51, 12077–12080.
12. Y.-S. Huang, H.-L. Xu, W.-J. Tian, Z.-S. Li, S. Escayola, M. Solà, A. Muñoz-Castro and Z.-M. Sun, J. Am. Chem. Soc., 2025, 147, 9407–9414.
13. N. M. Tam, H. T. Pham and M. T. Nguyen, Chem. Phys. Lett., 2014, 608, 255–263.
14. J. R. Leto, F. A. Cotton and J. S. Waugh, Nature, 1957, 180, 978–979.
15. S. Dorđević, J. Poater, M. Solà and S. Radenković, Chem. Sci., 2025, 16, 9920–9933.
16. S. Furukawa, M. Fujita, Y. Kanatomi, M. Minoura, M. Hatanaka, K. Morokuma, K. Ishimura and M. Saito, Commun. Chem., 2018, 1, 60.
17. B. Kiran, G. G. Kumar, M. T. Nguyen, A. K. Kandalam and P. Jena, Inorg. Chem., 2009, 48, 9965–9967.
18. M. Saunders, J. Comput. Chem., 2004, 25, 621–626.
19. F. Avaltroni and C. Corminboeuf, J. Comput. Chem., 2012, 33, 502–508.
20. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
21. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
22. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
23. T. Clark, J. Chandrasekhar, G. W. Spitznagel and P. V. R. Schleyer, J. Comput. Chem., 1983, 4, 294–301.
24. 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, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision C.01, Gaussian, Inc., Wallingford, CT, 2016.
25. P. W. Fowler, R. Zanasi, B. Cadioli and E. Steiner, Chem. Phys. Lett., 1996, 251, 132–140.
26. J. Jusélius, D. Sundholm and J. Gauss, J. Chem. Phys., 2004, 121, 3952–3963.
27. S. Lehtola, M. Dimitrova, H. Fliegl and D. Sundholm, J. Chem. Theory Comput., 2021, 17, 1457–1468.
28. G. Monaco, F. F. Summa and R. Zanasi, J. Chem. Inf. Model., 2021, 61, 270–283.
29. D. Y. Zubarev and A. I. Boldyrev, Phys. Chem. Chem. Phys., 2008, 10, 5207–5217.
30. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
31. T. Lu, J. Chem. Phys., 2024, 161, 082503.
32. L. Rincon, R. Almeida, J. E. Alvarellos, D. Garcia-Aldea, A. Hasmy and C. Gonzalez, J. Chem. Soc., Dalton Trans., 2009, 3328–3333.
33. A. Savin, B. Silvi, F. Coionna and F. Colonna, Can. J. Chem., 1996, 74, 1088–1096.
34. M. Kohout, Dgrid-5.0, 2017.
35. T. D. Goddard, C. C. Huang, E. C. Meng, E. F. Pettersen, G. S. Couch, J. H. Morris and T. E. Ferrin, Protein Sci., 2018, 27, 14–25.
36. R. C. Sabando, C. Riplinger, F. Wennmohs, F. Neese and G. Bistoni, J. Chem. Theory Comput., 2025, 21, 7920–7934.
37. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15, e70019.
38. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
39. D. Rappoport and F. Furche, J. Chem. Phys., 2010, 133, 134105.
40. T. B. Tai, P. Kadłubański, S. Roszak, D. Majumdar, J. Leszczynski and M. T. Nguyen, ChemPhysChem, 2011, 12, 2948–2958.
41. L. Van Duong, N. N. Tri, N. P. Hung and M. T. Nguyen, J. Phys. Chem. A, 2022, 126, 3101–3109.
42. X. Dong, S. Jalife, A. Vásquez-Espinal, J. Barroso, M. Orozco-Ic, E. Ravell, J. L. Cabellos, W. Y. Liang, Z. H. Cui and G. Merino, Nanoscale, 2019, 11, 2143–2147.
43. X. Dong, S. Jalife, A. Vásquez-Espinal, E. Ravell, S. Pan, J. L. Cabellos, W. Y. Liang, Z. H. Cui and G. Merino, Angew. Chem., Int. Ed., 2018, 57, 4627–4631.
44. B. Feng, J. Zhang, Q. Zhong, W. Li, S. Li, H. Li, P. Cheng, S. Meng, L. Chen and K. Wu, Nat. Chem., 2016, 8, 563–568.
45. Q. Li, E. B. Aklile, A. Tsui and M. C. Hersam, Nat. Chem., 2025, 17, 642–652.
46. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., 2009, 36, 199–207.
47. H. Fliegl, D. Sundholm, S. Taubert, J. Jusélius and W. Klopper, J. Phys. Chem. A, 2009, 113, 8668–8676.
48. J. C. Santos, W. Tiznado, R. Contreras and P. Fuentealba, J. Chem. Phys., 2004, 120, 1670–1673.

---