Triple-chain boron clusters: lengthening via phosphorus doping and enhancing stability through {P} by {CH} substitution

Abstract

This theoretical study presents novel insights into the doping of boron clusters with an increasing number of dopant atoms, ranging from 1 to 4, that preserve the integrity of the original boron framework. The triple-chain forms of clusters B₁₀ and B₁₆ remain unchanged upon sequential addition of P atoms, showcasing a perfect isolobal substitution of {P} with {CH}. Similarities in the number of delocalized electrons are observed between pure and doped boron clusters, alongside the subsequent substitution of {P} with {CH}. All triple-chain structures exhibit high thermodynamic stability, having low vertical attachment energies and high vertical ionization energies. The lowest-lying isomer of B₁₆P₄ has thus a triple-chain shape instead of a tetrahedral T_d form as previously reported. While AdNDP analysis confirms the number of globally delocalized electrons, it differs much from a previous interpretation. The magnetic ring current maps support the double aromaticity of triple-chain structures. Electron counting rules established for triple-chain structures are verified. The particle-in-a-rectangular-box model elucidates the relationship between the structure size and electron configuration and aids in understanding the transition from antiaromatic to aromatic configurations. The self-locking phenomenon is crucial for adhering to the triple-chain model and satisfying electron configuration requirements.

---

1. Introduction

The field of nanocluster science has continuously developed over the past few decades, and the growth mechanisms of the clusters of several elements have gradually emerged. Of the atomic clusters investigated, those formed by the boron atom exhibit the most diverse and non-conventional structural motifs and unpredictable geometry transformations following a small change in the size, net charge, and doping.¹ In a series of small pure clusters, going from B₃ to B₁₉, they are either planar or quasi-planar (2D)^2,3 with a fullerene-like intercalation of the size of B₁₄.⁴ The size of B₂₀ is renowned for its double ring (DR) tubular structure;⁵ after which the B_n class returns to the quasi-planar B₂₁ class.⁶ At larger sizes, there is still an alternating occurrence of the global minima between 2D and 3D structures.^7–12 Due to such a geometric diversity of boron clusters, various models of electronic structures have been proposed to rationalize their thermodynamic stability, and in particular their aromatic nature which in part governs their stability. Boron clusters appear to challenge the classical models of chemical bonding, requiring the introduction of numerous new electron count models that also necessitate a redefinition of the aromaticity concept.^13,14

For small size clusters, the classical Hückel rule proves effective for most planar and quasi-planar structures.^15,16 However, for circular structures such as the anion B₁₉⁻,^17,18 or the dianion B₂₀²⁻,¹⁹ the Hückel model becomes inadequate; in fact, it failed to rationalize their aromatic character. Instead, the disk aromaticity model^19–21 has proved to be more appropriate for the circular disk and even generalizable for a larger class of planar and bowl structures.

Instead of using the 2(N + 1)² electron count of the Hirsch rule²² or the phenomenological shell model (PSM),²³ the hollow cylinder model^24,25 provides us with a better understanding of the electronic configuration in tubular structures, especially when the tubes are formed from multiple atomic rings, making the electron picture more evident.²¹

The ribbon or double-chain motif^26–29 has emerged as a distinctive form due to its two boron chains, but because of impurities, often in a negatively charged state, it can now be extended up to a length of 17.0 Å.²⁹ The ribbon model^30,31 has subsequently been used to rationalize the aromaticity of ribbon-like structures without necessarily requiring cyclic electron delocalization¹³ in the structure. Surpassing mass spectrometry experiments for investigating boron ribbons in an anionic state,²⁷ the boron ribbon cluster B₂C₄H₄ was synthesized in a neutral state upon substitution of H atoms with appropriate ligands.³²

The structures formed from three boron chains³³ are less commonly found, and the electron counting rule of [σ²⁽ⁿ⁺¹⁾π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] was established for rationalizing their aromaticity.³⁴Scheme 1 presents two dianionic structures, B₁₀²⁻ and B₁₆²⁻, which represent the triple-chain motif best. The electron configuration of B₁₆²⁻ satisfies the [σ²⁽ⁿ⁺¹⁾π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] counting rule with n = 2. Scheme 1 also illustrates the global minima of the anionic and neutral B₁₀^−/0 and B₁₆^−/0,^2,35–37 whose shapes involve a slight geometric deformation from a perfect planarity of the corresponding dianions B₁₀²⁻ and B₁₆²⁻, respectively. In this context, a question of interest is whether it is possible to extend the size of boron triple-chain clusters by means of doping. Following multiple attempts to introduce various impurities to the latter, we find that the phosphorus element becomes both a neutralizing impurity and an extending agent. As a matter of fact, the incorporation of P atoms into the above dianionic species not only neutralizes the dianions but also extends the triple-chain form of the resulting B₁₀P₄ and B₁₆P₄ systems. A similar picture for the evolution of the structural motif of the B_xP_y clusters remains valid in the isovalent analogues in which each {P} atom is replaced by a {CH} unit (both entities having 5 valence electrons each). The substitution of {P} with {CH} ultimately leads to carborane structures that offer more practical opportunities for synthesis.

Scheme 1 Global minima of boron clusters B₁₀^2−/−/0 and B₁₆^2−/−/0 and the electronic configuration of the planar dianion B₁₆²⁻ including six σ and ten π orbitals.

The account presented in the next sections aims to contribute to a demonstration of the validity of the [σ²⁽ⁿ⁺¹⁾π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] electron counting rule for this class of triple-chain compounds. Our calculated results point out that the inherent electronic distribution, and thereby the corresponding aromatic character, favours the planar triple-chain form that emerges to be thermodynamically more stable than the more symmetrical tetrahedron or other three-dimensional isomers. This is particularly related to a recent theoretical study which reported that the tetrahedron is the lowest-lying isomer of the B₁₆P₄ cluster system.³⁸ In contrast, our results show that the most stable B₁₆P₄ isomer exhibits a perfect planar triple-chain geometry.

2. Computational methods

Isomers of each of the clusters considered in this theoretical study were located following full geometry optimizations of over 1000 different initial geometries. These initial geometries were generated using a stochastic search algorithm³⁹ and a genetic algorithm (GA)⁴⁰ that was implemented by our group.⁸ Geometries of the initial guess structures were optimized using density functional theory (DFT) with the PBE⁴¹ functional, which has recently been benchmarked and found to be suitable for examining the relative energies of clusters containing multiple B atoms,⁴² and the 6-31G(d)^43,44 basis set without subsequent vibrational frequency calculations. The set of isomers with relative energies within 3 eV with respect to the most stable structure is referred to as the seed structure ensemble. Each isomer in the seed structure ensemble was used to generate a series of input for optimization via both the stochastic search algorithm³⁹ and the genetic algorithm (GA).⁴⁰ After optimizations, the new set of output generated was evaluated to determine whether they belong to the seed structure ensemble. This iterative process continued until no new seed structures were identified.

The optimized structures whose total energies were within ∼3 eV from the most stable geometry (the final seed structure ensemble) were then reoptimized using the same PBE functional but with a larger 6-311+G(d)^43,44 basis set, along with vibrational frequency calculations in order to ensure the true local minimum character. Zero-point vibrational energy corrections were obtained from the harmonic frequencies without scaling. To identify the global minimum of each system, the electronic energies of a set of ten lowest-lying isomers were subsequently computed using single-point coupled-cluster theory CCSD(T)/6-311+G(d) calculations making use of the PBE/6-311+G(d) optimized geometries. Unrestricted Hartree–Fock (UHF) formalism was used for open-shell systems. All electronic structure calculations were performed using the Gaussian 16 software suite.⁴⁵

The adaptive natural density partitioning (AdNDP) algorithm,⁴⁶ as implemented in the Multiwfn software,⁴⁷ was used to probe the number of globally delocalized electrons. The chemical bonding and aromatic character were assessed through the electron localization function (ELF) maps^48,49 using the Dgrid 5.0 package.⁵⁰

Topological analysis of the electron localization function (ELF) map indicates that a structure with a high bifurcation value on the ELF isosurface displays aromatic characteristics, while that with a low bifurcation value does not exhibit aromaticity. A planar molecule often has aromaticity based on the behavior of its π and σ electrons. Santos et al. demonstrated that ELF_π and ELF_σ maps can help classify the aromaticity of these electron types.⁵¹ In certain specific cases, σ electrons can distinctly be categorized into localized σ electrons and delocalized σ electrons. At this point, the separation of ELF_σ into ELF_σl (for localized σ electrons) and ELF_σd (for delocalized σ electrons) further clarifies the impact of each electron type on the properties and stability of the structure considered.^30,52

Maps of the magnetic current density induced by a perpendicular magnetic field to the planar structure were calculated and visualized using the SYSMOIC⁵³ program; current flows were utilized to indicate the diatropic (clockwise) or paratropic (counterclockwise) nature of the ring current.

3. Results and discussion

3.1. B_xP_y systems

As stated above, we use phosphorus atoms as dopants in two boron systems displayed in Scheme 1. Four phosphorus atoms are thus successively added into the B₁₀ and B₁₆ skeletons. The lowest-lying isomers of the resulting B₁₀P_y and B₁₆P_y with y = 1, 2, 3 and 4 are presented in Fig. 1 and 2. The relative energy ordering of isomers at 0 K established by (U)CCSD(T) calculations and corrected for zero-point energies (ZPEs) are marked by the labels having increasing numbers of 1, 2, 3, etc…

Fig. 1 Optimized geometries of the lowest-lying B₁₀P_y isomers with y = 1, 2, 3, and 4. Relative energies, given in kcal mol⁻¹, are obtained from (U)CCSD(T)/6-311+G(d) + ZPE computations; PBE/6-311+G(d) + ZPE values are given in parentheses.

Fig. 2 Optimized geometries of the lowest-lying B₁₆P_y isomers with y = 1, 2, 3, and 4. Relative energies, given in kcal mol⁻¹, are obtained from (U)CCSD(T)/6-311+G(d) + ZPE computations; PBE/6-311+G(d) + ZPE values are given in parentheses.

Fig. 1 and 2 illustrate two prominent trends when phosphorus atoms are added into the boron triple-chain structures. Firstly, P atoms are preferentially attached to positions 1, 2, 3 and 4 of four B–B bonds at both ends of the structures, as marked in Scheme 2. The next preferred positions are the B–B bonds denoted as 5, 6, 7 and 8 on the outermost of the two outer B rows.

Scheme 2 Identification of the preferred binding sites of P atoms into the B–B bonds of the B_10+6n triple-chain framework.

The second trend is a rapid flattening of the triple-chain structure following doping, where, upon the second P atom addition, the triple-chain form already becomes perfectly planar, and such a perfect planarity is maintained following doping with the 3^rd and 4^th P atoms. This is particularly intriguing as it basically differs from a typical observation on boron cluster doping which tends to favour 3D isomers. A recent study conducted by Jin et al.⁵⁴ reported that the doping of an aluminium atom into the B₁₀ size leads to the most stable 3D isomer in which the Al atom is positioned on the surface of a triple-chain structure, while subsequent lower-lying isomers of AlB₁₀ do not feature any triple-chain framework. Similarly, stable isomers of AlB₁₆ also do not have any triple-chain skeleton.⁵⁴ A previous study by Mai and co-workers⁵⁵ indicated that the silicon counterpart of B₁₀P-3 (Fig. 1) is the global minimum of the SiB₁₀ cluster in which the Si atom is attached to one long edge of the triple-chain structure, while the triple-chain framework is not even present in higher-energy isomers. A subsequent study by Mai and co-workers⁵⁶ also demonstrated that SiB₁₆ has a second lowest-lying isomer in which the Si atom lies in the plane of a triple-chain B₁₆, while the global minimum and other isomers lack a triple-chain skeleton.

With our current knowledge, no previous observation was reported on the doping effect of an element, excluding H, halogens and noble gases, into a small-sized boron cluster whose resulting doped structure retains the original pure boron framework by up to four dopant atoms. The 3D structure that appeared in B₁₀P₄ is B₁₀P₄-4 which is 14.8 kcal mol⁻¹ higher in energy than the global minimum B₁₀P₄-1. In the case of B₁₆P₄, the 3D structure bearing the form of a T_d tetrahedron B₁₆P₄-2 is only 5.6 kcal mol⁻¹ higher than B₁₆P₄-1.

The T_d tetrahedron P₄ is well known for its superior thermodynamic stability,⁵⁷ which is believed to result from its electron shell satisfying the doubly spherically aromatic count of 2(N + 1)² with 18σ and 2π electrons.⁵⁸ Thus, insertion of a B₁₆ unit into the tetrahedral P₄ forming the B₁₆P₄-2 isomer is highly noteworthy. The existence of both B₁₆P₄-1 and B₁₆P₄-2 has been recently reported by Gao et al.³⁸ In their study, DFT computations using the PBE0 functional suggested that B₁₆P₄-2 (T_d) is more stable than B₁₆P₄-1, despite the results from TPSSh, or more importantly, CCSD(T) indicates a reversed energy ordering.

In the present study, we find agreement in the energy ordering between the PBE functional and CCSD(T) results, indicating that in the neutral state B₁₆P₄-1 is more stable than B₁₆P₄-2. Based on PBE functional calculations, B₁₆P₄-1 is more stable than B₁₆P₄-2 by 27 kcal mol⁻¹ in the anionic state. Consequently, it is more likely that the B₁₆P₄-1 can experimentally be detected.

Our present result disagrees with a previous report by Gao et al.,³⁸ which concluded that the T_d structure of B₁₆P₄-2 is the most stable isomer. These authors³⁸ also relied on the adaptive natural density partitioning (AdNDP) analysis to explain the stability of B₁₆P₄-2. In addition, the low-lying B₁₆P₄-3 and B₁₆P₄-4 isomers were missed in their study. To have another look at these conflicting results, we replicate the analysis made in ref. 38, the details of which are given in the ESI file.†

3.2. B_xP_yC_zH_z systems

The concept of “isolobal analogy”, meaning “similar lobes”, developed by Hoffmann⁵⁹ proved to be helpful in designing novel molecules based on stable molecular structures, and numerous molecules have been discovered using this approach.^60,61 Considering the fact that a {P} atom is not only isovalent but also isolobal to a {CH} fragment, we now carry out a sequential substitution of {P} atoms by {CH} units, and the corresponding lowest-lying isomers are presented in Fig. 3 and 4. Furthermore, in the anionic state, the triple chain of B_xP_yC_zH_z retains the global minima, as illustrated in Fig. S2 and S3 (ESI file†).

Fig. 3 Optimized geometries of the lowest-lying B₁₀P_yC_zH_z isomers with y = 3, 2, 1 and 0, and z = 4 − y. Relative energies given in kcal mol⁻¹ are obtained from CCSD(T)/6-311+G(d) + ZPE computations; PBE/6-311+G(d) + ZPE values are given in parentheses.

Fig. 4 Optimized geometries of the lowest-lying B₁₆P_yC_zH_z isomers with y = 3, 2, 1 and 0, and z = 4 − y. Relative energies, given in kcal mol⁻¹, are obtained from PBE/6-311+G(d) + ZPE computations.

Calculated results indicate the emergence of several stable isomers with a close relative energy (RE) as a consequence of the multiple permutations of positions between {P} and {CH}. The 1-4-isomer is the isomer with {P} or {CH} bonds at the 1, 2, 3 and 4 positions as shown in Scheme 2, and the 5-(2-4)-isomer is that with {P} or {CH} bonds at the 2, 3, 4 and 5 positions. Even more interestingly, we observe a gradual increase in the RE of the most stable isomer (1-4-isomer) as compared to the 5-(2-4)-isomer. This RE increases from 3.0 to 7.9 kcal mol⁻¹ going from B₁₀P₄ to B₁₀C₄H₄. From B₁₆P₄ to B₁₆C₄H₄, the RE increases from 6.8 to 9.6 kcal mol⁻¹, and eventually the triple-chain B₁₆C₄H₄-1 exhibits a higher likelihood for experimental synthesis with a minimal occurrence of other isomers. The recent synthesis of the ribbon-like carborane B₄C₂R₄ ³² constitutes a stimulation for a promising generation of these triple-chain carboranes B_xC₄H₄, which in turn could facilitate the synthesis of B_xP_yC_zH_z structures in a reversed manner.

The isomers derived from the T_d structure of B₁₆P_xC_yH_y presented in Fig. 4 emphasize a rapid increase in the REs between both triple-chain and T_d forms when {P} is sequentially substituted by {CH}. This is analogous to the transition from the 3D to the 2D form observed in the tetrahedral P₄ when {P} is successively replaced by {CH} forming C₄H₄.^62,63

3.3. Stability trend

The vertical attachment energy (VAE) and the vertical ionization energy (VIE) are defined as:

VAE = E(X⁰) − E(X⁻) (1)

VIE = E(X⁰) − E(X⁺) (2)

where E(X^−/0/+) is the energy of cluster X in anionic, neutral and cationic states, respectively, at the optimized neutral geometry. The chemical hardness⁶⁴ (η), electronegativity⁶⁵ (χ) and electrophilicity⁶⁶ (ω) are then calculated using formulas (3), (4) and (5):where the chemical potential (μ) is the opposite of electronegativity (μ = −χ).

The frontier orbital HOMO–LUMO gaps (HLGs) of B_xP_yC_zH_z are presented in Fig. 5, whereas the VAE, VIE, η, χ and ω values are plotted in the graphs shown in Fig. 6, and given in Tables S1 and S2 of the ESI file.†

Fig. 5 Energy diagram of triple-chain B_xP_yC_zH_z. HOMO–LUMO gaps (eV) are given in grey areas.

Fig. 6 VAE, VIE, η, χ and ω values of B_xP_yC_zH_z clusters in (a) x = 10 and (b) x = 16.

The HLG and η values point toward increasing structural stability as their values increase, whereas ω indicates that stability increases as its values fall. A high VIE suggests that the structure tightly holds its electrons, whereas a low VAE indicates a stabilized neutral state. Therefore, a high VIE and low VAE values show the stability of a structure in its neutral state. Overall, while the HLG, VAE and η values of B₁₆P_yC_zH_z are lower than those of B₁₀P_yC_zH_z, VIE and ω of B₁₆P_yC_zH_z are higher than that of B₁₀P_yC_zH_z; B₁₀P_yC_zH_z thus have greater stabilization as compared to B₁₆P_yC_zH_z.

Initially, when P is doped into B_x, the values of HLG, VAE, VIP, η and ω vary according to parity; closed-shell systems are more kinetically stable than their open-shell counterparts. Considering only the closed-shell systems, their stability decreases upon P addition: B_x is more stabilized than B_xP₂, and B_xP₂ is more stable than B_xP₄. Subsequently, the replacement of {P} with {CH} leads to the opposite effect: the more {P} is replaced by {CH}, the more stable the structure becomes.

The electronegativity (χ) of a molecule refers to its ability to attract electrons towards itself. The electronegativity of structures considered here exhibits less variability as compared to the aforementioned parameters with the lowest value being 5.0 eV for B₁₆C₄H₄ and the highest value being 5.7 eV for B₁₀.

3.4. Chemical bonding

The AdNDP method is often employed with the aim of verifying the nature of the MOs in the species considered. The set of σ bonds of the B₁₆²⁻ dianion, which was assigned from AdNDP analysis carried out by Sergeeva et al.,³⁷ suggests that B₁₆²⁻ has eight σ delocalized electrons leading to its σ antiaromatic character, according to the classical (4n) electron count. Our subsequent study,³⁴ which was a thorough analysis of the σ-out-of-plane contributions to the molecular magnetizability, revealed that B₁₆²⁻ rather exhibits σ aromaticity with the sum of magnetic components . Our present AdNDP analysis is summarized in Fig. 7 whereas the corresponding MOs are shown in Scheme 1. The ring current maps for σ electrons of B₁₆²⁻ displayed in Fig. S8 (ESI file†) also confirm the σ aromatic character of this dianion.

Fig. 7 Chemical bonding pattern of the B₁₆²⁻ dianion based on AdNDP analysis. Occupation numbers (ONs) are indicated. (a) is the set of σ bonds obtained by replicating Sergeeva's results (ref. 37) excluding bond (7); (b) is the set of σ bonds obtained by the method employed in this study, and (c) represents a set of π bonds, applicable to both (a) and (b).

In view of such a significant discrepancy, we conducted the AdNDP analysis again, and our recent results are presented in Fig. 7a, corresponding to the σ bonding framework of B₁₆²⁻. The set of eighteen σ bonds from (1) to (6) are well reproduced, similar to those reported by Sergeeva et al.,³⁷ despite a slight discrepancy in the occupation numbers (ONs) that are rather not significant due to the utilization of different theoretical methods. It is noteworthy that after collecting eighteen σ bonds (cf. (1–6) in Fig. 7a) and five π bonds (Fig. 7c), a significant amount of residual valence electrons remains sufficient to be considered as a bond. Indeed, a σ 16c–2e bond (7) is obtained from this set of residual valence electrons with an ON value of 1.80 |e|, demonstrating its delocalized nature. This points out that the analysis reported in ref. 37 is not appropriate, thereby prompting us to repeat it by setting the occupation threshold to 1.9 |e| to first identify the bonds with ON values being larger than 1.9 |e|.

The σ bond framework as depicted in Fig. 7b reveals that the four 2c–2e σ bonds (3) turn out to be four 3c–3e σ bonds (10) with ON values now exceeding 1.9 electrons, the six σ 3c–2e bonds (4–5) turn out to be only three σ 4c–2e σ 3c–2e bonds (12–13), and finally, the set of remaining residual valence electrons comprises three globally delocalized σ bonds with shapes analogous to those of the σ delocalized MOs indicated in Scheme 1 (or Fig. S4a in the ESI file†).

The π bond framework also contains five globally delocalized π bonds whose shapes are akin to the π MOs as indicated in Scheme 1. Overall, the B₁₆²⁻ dianion truly comprises 3 globally delocalized σ MOs and 5 π MOs as illustrated in Scheme 1. The remaining σ electrons include localized (2c–2e) electrons and restricted delocalized (nc–2e) electrons with n = 3 or 4. Such a chemical bonding pattern, which is similar to those of the B₂₂²⁻, B₁₀P₄, B₁₀C₄H₄, B₁₆P₄ and B₁₆C₄H₄ systems, is presented in Fig. S5, S6 and S7 (ESI file†), and demonstrates the assignment of the MOs pictured in Fig. 8 and Fig. S4 (ESI file†) as globally delocalized MOs are entirely justified.

Fig. 8 The σ delocalized and π MOs of (a) B₂₂²⁻ and (b) B₁₆P₄. H and L stand for the HOMO and LUMO, respectively.

The electron counting rule of [σ²⁽ⁿ⁺¹⁾π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] which has been proposed in our previous studies,³⁴ appears to be suitable for rationalizing the aromatic character of the triple-chain B₁₀²⁻ and B₁₆²⁻, and the antiaromaticity of the perfectly planar triple-chain B₂₂²⁻.

It should be reiterated again that this electron counting rule stems from consideration of the π electrons of the triple-chain structure as freely moving within a rectangular potential well with infinitely high potential at the boundaries. The corresponding Schrödinger equation for a particle moving in a rectangular box is a well-known basic concept found in many textbooks of quantum chemistry. Its eigenvalues are presented in eqn (6) that express the energy levels:

in which a and b are the length and the width of the rectangular box, respectively, and p = 1, 2, 3, … and q = 1, 2, 3, … are the length and width quantum numbers, respectively. While the wavefunction ψ^π_(p,1) with q = 1 belongs to the π₁ group, the wave function ψ^π_(p,2) with q = 2 belongs to the π₂ group.

The σ and π MOs of both B₁₆²⁻ and B₂₂²⁻ as well as those of the doped B₁₀P₄ and B₁₆P₄ are shown in Fig. 8 and Fig. S4 (ESI file†). Here, the electron configuration [σ⁶π₁⁶π₂⁴] of B₁₀P₄ coincides with that of B₁₆²⁻ and indicates that B₁₀P₄ achieves triple-chain aromaticity.

The B₂₂²⁻ dianion exhibits triple-chain antiaromatic characteristics due to the electron configuration of [σ⁸π₁¹⁰π₂⁴]. Interestingly, there is a permutation between the two wavefunctions ψ^π_(5,1) and ψ^π_(3,2), where ψ^π_(5,1) is the HOMO in B₂₂²⁻ and becomes the LUMO in B₁₆P₄, and ψ^π_(3,2) is the LUMO in B₂₂²⁻ and becomes the HOMO−1 in B₁₆P₄. This leads to fulfilment of the electron counting rule of [σ⁸π₁⁸π₂⁶] for B₁₆P₄. According to this triple-chain model, both B₁₀P₄ and B₁₆P₄ behave as triple-chain aromatic species.

Each of the B₂₂²⁻ and B₁₆P₄ species exhibits 14 π electrons and thus appears to be aromatic according to the classical (4n + 2) rule with n = 3. As pointed out in our previous study³⁴ and as also demonstrated by the π ring current maps illustrated in Fig. 9, while B₂₂²⁻ turns out to be a π antiaromatic species, B₁₆P₄ is a π aromatic one. Thus, the classical (4n + 2) electron counting rule does not account for the aromaticity of a triple-chain structure, whereas the [π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] electron count represents the more suitable rule related to the triple-chain aromaticity, albeit its representation is somewhat more complex.

Fig. 9 The σ and π ring current maps of B₂₂²⁻ and B₁₆P₄ triple-chain structures.

The (4n + 2) rule tends to lead us to associate these electrons with one non-degenerate level and n doubly degenerate levels of energy. However, many related structures are characterized with an electron count satisfying the (4n + 2) rule but they do not have any degenerate energy level, such as the oligoacenes. Several recently discovered structures exhibit aromatic properties and are explained by having electron spectra similar to that of naphthalene (C₁₀H₈)^37,67 or anthracene (C₁₄H₁₀).⁶⁸ Such a comparison is only partially correct, because while the π electron configuration of boron clusters can be compared to those of hydrocarbons on the basis of the similarity between the π-MO shapes and the number of π electrons, the boron clusters have in addition a set of σ delocalized electrons that do not exist in hydrocarbons. Accordingly, it appears better to avoid an explanation for the aromaticity of boron clusters using the classical (4n + 2) electron counting rule. Indeed, the electron configurations of oligoacenes including naphthalene, anthracene and tetracene that are presented in Fig. S11 (ESI file†) demonstrate well that they follow the [π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] rule. In summary, the stability of boron clusters often needs to be understood considering both sets of π and σ characteristics, whereas hydrocarbons only involve their π electrons.

Two additional issues need some further attention, namely, why triple-chain structures need to satisfy the [π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] rule, and how a cluster structure can satisfy this rule. Answers for both issues have been presented in our previous study,³⁴ but here, taking the doping effect into account, we would again demonstrate that the above rationalization is entirely valid.

The initial answer can be found through the relationship between the electron localization function (ELF) maps and bond lengths of structures as depicted in Fig. 10 and S10 (ESI file†). Both ELF_σd and ELF_π maps represent the ELFs originating from both sets of globally σ and π delocalized electrons. In these figures, the structures are highlighted with purple and yellow rhombuses. Purple rhombuses correspond to locations where multiple interactions of ELF_σd and ELF_π basins are present, resulting in shorter edge lengths and diagonals as compared to the edges and diagonals of the yellow rhombuses.

Fig. 10 The ELF_σd, ELF_π, and the structures of (a) B₂₂²⁻ and (b) B₁₆P₄. Bond lengths are given in angstrom.

The first observation is that in ELF_σl (cf. Fig. S9, ESI file†), which is similar to what has been observed in ribbon-like,³⁰ double ring,⁵² or triple ring²⁴ structures, the basins of attraction are only observed for bonds within the same row (sr-bond), not connecting rows (cr-bond), except for the bonds at the ends of each structure. Therefore, the common resulting characteristic is that the lengths of sr-bonds are shorter than those of cr-bonds. The ELF_π basins also indicate that they enhance more the bonding of sr-bonds than that of the cr-bonds. Within the same row, the bonds where ELF_π basins are present tend to have shorter distances than those of the same type without ELF_π basins.

The B2–B4 bond length (cf. Fig. S10 ESI file† for atom labels) corresponds to the shortest B–B bond that amounts to 1.54, 1.61 and 1.60 Å in B₁₆²⁻, B₁₀P₄ and B₁₀C₄H₄, respectively. The B4–B6 (B2–B9) bonds are the smallest diagonals within the yellow rhombus adjacent to the central purple rhombus. The deviation between B4–B6 and B2–B4 in B₁₆²⁻, B₁₀P₄ and B₁₀C₄H₄ amounts to 0.09, 0.17 and 0.24 Å, respectively. This significant difference is analogous to the presence of the ELF_π basin in the B2–B4 bond while such a basin is absent in the B4–B6 bond. In B₁₆²⁻, the B6–B16 bond distance is 1.6 Å, and it is also shorter than the neighbouring B4–B6 bond for the same reason. Thus, ELF_σd, in turn, tends to show an enhancement of the cr-bonds.

An interesting observation from these results is that the rhombuses with small diagonals that are enhanced continue to be reinforced by ELF_σd, while the rhombuses with small diagonals that are not enhanced are not reinforced by ELF_σd. We present herewith the clearest observation of this trend. The B6–B7 bond in B₁₆²⁻ is still influenced by ELF_σd, resulting in a bond length of 1.80 Å. However, in B₁₀P₄ it is no longer affected by ELF_σd, and this bond length is then increased to 1.85 Å. In B₁₀C₄H₄ the latter distance further increases to 1.91 Å. The B6–B7 bond in B₁₀P₄ corresponds to the B12–B13 bond in B₁₆P₄ and B₁₀C₄H₄ whose respective bond lengths amount to 1.85 and 1.90 Å.

There is an alternating arrangement of domains, coloured as purple and yellow rhombuses, in the triple-chain structure when it achieves aromaticity (cf.Fig. 10). Indeed, such alternation does not exist in B₂₂²⁻ and therefore the B₂₂²⁻ dianion is considered as a π antiaromatic species.

Although both B₁₆P₄ and B₂₂²⁻ have the same number of σ and π electrons, only B₁₆P₄ is aromatic. The enhanced and non-enhanced rhombuses observed in B₁₆P₄ indicate that the triple-chain aromaticity requires its electrons to be distributed in such a way that there is an alternating arrangement of rhombuses in the structure, as described above. This alternation also appears in triple-chain hydrocarbons as previously presented for the ribbon clusters.³⁴

The alternation between enhanced and non-enhanced regions in triple-chain structures³⁴ or ribbon structures³⁰ challenges the common notion of bond length equalization⁶⁹ for an aromatic structure (such as in the case of benzene). This challenge is further reinforced when the B₂₂²⁻ structure, which exhibits even more bond length equalization, is now identified as a π antiaromatic species.

In the ELF_σd and ELF_π maps (cf.Fig. 10), large anomalous basins have appeared that are highlighted in green. Due to the insignificant difference between the bonds affected by these basins in B₁₆²⁻ as compared to the adjacent bonds, the role of these basins was not mentioned in our previous study.³⁴ Because of the strong influence of impurities, the roles of these basins become more apparent. The significantly shorter bond length of P12–B7 (C12–B7) as compared to P12–B6 (C12–B6) in B₁₀P₄ (B₁₀C₄H₄) or P18–B13 (C18–B13) compared to P18–B12 (C18–B12) in B₁₆P₄ (B₁₆C₄H₄) suggests a self-locking phenomenon akin to the ribbon model.³⁰ Indeed, these B–P (or B–C) bonds that are anomalously short with respect to their normal lengths, and are located at the four corners of the rectangular structure, indicate that the aromatic motif in an triple-chain model self-adjusts to the particle in a rectangular box (PIARB) model, which has infinite potential at the boundary, making it the most suitable electron configuration possible.

Although they have the same number of π electrons, the correspondence between the antiaromatic B₂₂²⁻ to the aromatic B₁₆P₄ clearly demonstrates the significance of the triple-chain model. In B₂₂²⁻, the ψ^π_(5,1) level is occupied and the ψ^π_(3,2) level is unoccupied, whereas in B₁₆P₄ the ψ^π_(3,2) level is occupied and the ψ^π_(4,1) level is unoccupied (cf.Fig. 8). That means, if the ψ^π_(p2,2) level is the molecular orbital with p2 being the highest quantum number for the π₂ group, then the energy level of ψ^π_(p2+2,2) must be higher than that of ψ^π_(p2,2), and the structure must have enough π electrons for the ψ^π_(p2+1,1) level to be fully occupied. If we set p1 = p2 + 2 and according to eqn (6), the above condition implies:

orwith the SR being the standard a/b ratio for which a/b must be smaller, in such a way that the energy of ψ^π_(p2+2,1) is higher than that of ψ^π_(p2,2).

In a recent study, a small and simple Python code snippet was written to determine the dimensions a and b of the rectangle in which the π electrons move freely,⁴² in which the root mean square error (RMSE) index was calculated for determining the accuracy of the model (eqn (8)):

with m being the number of occupied molecular orbitals; E^DFT_i and E^PIARB_i represent the energies of the MOs calculated from the DFT method and from the PIARB, respectively.

Applying this code, the results for the structures C₁₀H₈, B₁₆²⁻, B₁₀P₄ and B₁₀C₄H₄ for 10 π electron systems and C₁₄H₁₀, B₂₂²⁻, B₁₆P₄ and B₁₆C₄H₄ for 14 π electron systems are provided in Table 1. With the RMSE < 0.01 for all systems considered in Table 1, the π electrons appear to follow the PIARB model.

Table 1 Length a and width b of C₁₀H₈, B₁₆²⁻, B₁₀P₄, B₁₀C₄H₄, C₁₄H₁₀, B₂₂²⁻, B₁₆P₄ and B₁₆C₄H₄ and their distance ratio a/b, standard ratio SR, and the RMSE value of each structure

Species	a (Å)	b (Å)	a/b	SR	RMSE
C10H8	9.21	6.52	1.41	2.00	0.0052
B16^2−	10.61	5.75	1.85	2.00	0.0022
B10P4	10.27	6.90	1.49	2.00	0.0065
B10C4H4	10.34	6.97	1.48	2.00	0.0052
C14H10	11.98	6.57	1.82	2.31	0.0057
B22^2−	13.84	5.78	2.39	2.31	0.0021
B16P4	13.55	6.54	2.07	2.31	0.0056
B16C4H4	13.56	6.55	2.07	2.31	0.0044

---

The heteromolecular compounds with an RMSE being twice as large as that of the pure boron triple-chain clusters indicate the presence of the phenomenological shell model (PSM)²³ on the electron spectrum, but it remains non-significant. As the triple-chain size increases, the number of π electrons also increases, proportional to the a/b ratio and SR. However, the ratio a/b is found to increase more rapidly than the SR, leading a/b to approach the SR or even exceed the SR as in B₂₂²⁻. Substitution of a {B₆²⁻} unit with a {P₄} atom does not only decrease the length a but also increases the width b of the triple-chain form, thereby reducing the a/b ratio. In other words, there is no significant change in the variation of a and b following substitution of {P} with {CH}.

As the triple-chain size increases, the levels ψ^π_(p2,2) and ψ^π_(p2+2,1) in the oligoacenes approach each other more closely, meaning that the frontier orbital (HOMO–LUMO) energy gap (HLG) decreases, thereby increasing the multi-reference nature of the molecule and reducing its kinetic stability.^70–72 The perfect planarity of B₂₂X₄ with X = {P} or {CH} is also due to an electron configuration satisfying the [σ²⁽ⁿ⁺¹⁾π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] rule, but the HLG amounts to only 0.13 and 0.31 eV for B₂₂P₄ and B₂₂C₄H₄, respectively, leading to the appearance of imaginary frequencies for high symmetry structures. Therefore, these clusters are not further discussed in the present work. Furthermore, in boron-based clusters, σ electrons can also serve as the HOMO as seen in B₁₆P₄ or the LUMO as in the B₁₅⁻ anion,⁴² thus making the prediction of the size of the HLG as well as their aromatic character and stability even more challenging.

For a broader perspective on the effect of substitution of {P} with {CH}, the energy diagrams of the structures studied are presented in Fig. 5. At a rapid glance, substitution of {P} with {CH} leads to an increase in HLG, corresponding to enhanced kinetic stability as indicated above. The increase in HLG stems from two opposing processes resulting from the substitution, namely the downward shift of σ MOs and the upward shift of π MOs. This phenomenon is observed more distinctly in B₁₀P_yC_zH_z as compared to B₁₆P_yC_zH_z. In B₁₀P_yC_zH_z, this displacement leads to a rearrangement of the HOMO and HOMO−1 cross positions in B₁₀PC₃H₃. Therefore, under strict consideration, only B₁₀P₄, B₁₀P₃CH and B₁₀P₂C₂H₂ are deemed to be isolobal, while the entire set of B₁₀P_yC_zH_z structures are better considered as isovalent.

In B₁₆P_yC_zH_z, the ψ^σ₄ level becomes HOMO−1 in B₁₆PC₃H₃ and continues to shift downwards with increasing substitution of {P} with {CH}. Therefore, all B₁₆P_yC_zH_z isomers are considered as isolobal, except for B₁₆P₄. Both the HOMO and LUMO of B₁₆P_yC_zH_z, except for B₁₆P₄, are π MOs, and their energy levels increase. Consequently, the increase in HLG for B₁₆P_yC_zH_z is much slower as compared to that of B₁₀P_yC_zH_z. The increase in HLG for B₁₆P_yC_zH_z is also facilitated by a faster increase in the energy levels of LUMO (antibonding) as compared to HOMO (bonding).

With their π aromaticity, the triple-chain B_xP_yC_zH_z that are analogues of the aromatic oligoacenes and further enhanced by intrinsic σ aromaticity, along with a wide range of HLG spectra spanning from 1.1 to 2.2 eV, offer a diverse selection for experimental studies, with tailored applications according to specific purpose on triple-chain boron clusters. To emphasize further that boron triple-chain structures are not entirely analogous to acenes, the UV-Vis spectra of boron triple-chain compounds and oligoacenes displayed in Fig. S12 and S13,† with identical π electron configurations are presented in Fig. S10 (ESI file†). These electronic spectra demonstrate that triple-chain boron compounds exhibit the ability to absorb light across a broad range in the visible region, whereas oligoacenes mainly absorb light within a much narrower range and exhibit a pronounced blue-shift in the ultraviolet region.

4. Concluding remarks

The results obtained in the present theoretical study offer us a case of doped boron clusters in which the original boron framework remains largely unchanged as the number of phosphorus dopant atoms goes from 1 to 4. Indeed, both B₁₀ and B₁₆ host systems maintain their triple-chain form following sequential addition of P atoms. Furthermore, our calculated results indicate a perfect isovalent/isolobal substitution of {P} with {CH}. The similarity in the number of delocalized electrons formed between the triple-chain form of the pure boron B_n²⁻ and the doped B_n−6P₄ clusters and in those formed by the substitution of {P} with {CH} thereafter are intriguing. A remarkable effect of P substitution is neutralization of both dianionic hosts. In each series, the triple-row form remains by far the global minimum, ensuring a high likelihood of its experimental detection.

This study does not diminish the significance of AdNDP in the analysis of chemical bonding in molecules. However, the subjective nature of AdNDP can lead to user-dependent results. Relying solely on AdNDP analysis can result in a significant misunderstanding of the chemical bonding involved. Therefore, it is essential to employ additional methods to verify the AdNDP findings.

The number of globally delocalized σ and π electrons present in these clusters and the magnetic ring current maps contribute to the verification that B_xP_yC_zH_z triple-chain structures exhibit double aromaticity.

While the electron counting rule for the double-chain and ribbon species is [σ²⁽ⁿ⁺¹⁾π²ⁿ], the electron counting rule of [σ²⁽ⁿ⁺¹⁾π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] is established in the present study for the triple-chain form. The latter is first introduced for the dianions B₁₀²⁻ and B₁₆²⁻, and following doping, it is confirmed to be accurate for B₁₀P₄ and B₁₆P₄ through a similarity in electron configuration between the dianions and the doped derivatives. It can also be verified going from the configuration of [σ⁸π₁¹⁰π₂⁴] in B₂₂²⁻ to that of [σ⁸π₁⁸π₂⁶] in B₁₆P₄.

The model of a particle moving in a rectangular box proved to be effective in elucidating the relationship between the size of a structure and its electron configuration, thereby explaining the transition from an antiaromatic configuration of B₂₂²⁻ to an aromatic configuration of B₁₆P₄. The model of a particle moving in a rectangular box also helps us rationalize the decrease in the frontier orbital energy gap (HLG) of oligoacenes as their size increases.

The electron localization function (ELF) maps and their relationship with different types of bond lengths highlight the self-locking phenomenon as a crucial point for a cluster structure to adhere to the triple-chain model and the [σ²⁽ⁿ⁺¹⁾π₁²⁽ⁿ⁺¹⁾π₂²ⁿ] electron configuration to satisfy the regions of enhanced rhombuses and alternately non-enhanced rhombus areas.

All structures considered have vertical attachment energies (VAEs) smaller than 4 eV and vertical ionization energies (VIEs) larger than 7 eV. The substitution of a {P} atom with a {CH} unit indeed enhances the thermodynamic stability of the boron triple-chain form, which helps facilitate its laboratory synthesis and open avenues for a new class of materials with versatile electronic properties. The calculated UV-Vis spectrum also suggests its potential application as a semi-conductor.

Data availability

AdNDP analysis for the T_d tetrahedron B₁₆P₄-3 isomer and triple-chain structures; the σ delocalized and π-MOs of triple-chain structures and π-MOs of some oligoacenes; the σ and π ring current maps of triple-chain structures; ELF_σl, ELF_σd, and ELF_π of triple-chain structures; the predicted electronic absorption spectrum of triple-chain and some oligoacenes; and tables of the HOMO–LUMO gap (HLG), VIE, VAE, hardness (η), electronegativity (χ) and electrophilicity (ω) of triple-chain structures.

The quantum chemical program Gaussian 16 has been used. The reference to this program is given in the list of references. The datasets supporting this article are uploaded as part of the ESI file† that also includes the Cartesian coordinates of the optimized geometries and information on their electronic structures.

Conflicts of interest

The authors declare no conflicts of interest.

Acknowledgements

LVD and MTN are grateful to Van Lang University. This work was undertaken with the assistance of resources to AK from the National Computational Infrastructure (NCI), which is supported by the Australian Government.

References

1. J. Barroso, S. Pan and G. Merino, Chem. Soc. Rev., 2022, 51, 1098–1123.
2. T. B. Tai, D. J. Grant, M. T. Nguyen and D. A. Dixon, J. Phys. Chem. A, 2010, 114, 994–1007.
3. B. J. Nagare, S. Chavan and V. Bambole, Comput. Theor. Chem., 2018, 1125, 54–62.
4. L. Cheng, J. Chem. Phys., 2012, 136, 104301.
5. B. Kiran, S. Bulusu, H.-J. H.-J. Zhai, S. Yoo, X. C. Zeng and L.-S. L.-S. Wang, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 961–964.
6. R. Casillas, T. Baruah and R. R. Zope, Chem. Phys. Lett., 2013, 557, 15–18.
7. T. B. B. Tai and M. T. T. Nguyen, Phys. Chem. Chem. Phys., 2015, 17, 13672–13679.
8. H. T. Pham, L. Van Duong, B. Q. Pham and M. T. Nguyen, Chem. Phys. Lett., 2013, 577, 32–37.
9. X. Wu, L. Sai, S. Zhou, P. Zhou, M. Chen, M. Springborg and J. Zhao, Phys. Chem. Chem. Phys., 2020, 22, 12959–12966.
10. J. Lv, Y. Wang, L. Zhu and Y. Ma, Nanoscale, 2014, 6, 11692–11696.
11. T. B. Tai and M. T. Nguyen, Nanoscale, 2015, 7, 3316–3317.
12. A. B. Rahane and V. Kumar, Nanoscale, 2015, 7, 4055–4062.
13. G. Merino, M. Solà, I. Fernandez, C. Foroutan-Nejad, P. Lazzeretti, G. Frenking, H. L. Anderson, D. Sundholm, F. Cossio, M. A. Petrukhina, J. Wu, J. Wu and A. Restrepo, Chem. Sci., 2023, 21, 1–9.
14. H. Ottosson, Chem. Sci., 2023, 14, 5542–5544.
15. A. E. Kuznetsov and A. I. Boldyrev, Struct. Chem., 2002, 13, 141–148.
16. Q. S. Li and H. W. Jin, J. Phys. Chem. A, 2002, 106, 7042–7047.
17. W. Huang, A. P. Sergeeva, H.-J. Zhai, B. B. Averkiev, L.-S. Wang and A. I. Boldyrev, Nat. Chem., 2010, 2, 202–206.
18. J. O. C. Jiménez-Halla, R. Islas, T. Heine and G. Merino, Angew. Chem., Int. Ed., 2010, 49, 5668–5671.
19. T. B. Tai, R. W. A. Havenith, J. L. Teunissen, A. R. Dok, S. D. Hallaert, M. T. Nguyen and A. Ceulemans, Inorg. Chem., 2013, 52, 10595–10600.
20. T. B. Tai, L. Van Duong, H. T. Pham, D. T. T. Mai and M. T. Nguyen, Chem. Commun., 2014, 50, 1558–1560.
21. L. Van Duong and M. T. Nguyen, Reference Module in Chemistry, Molecular Sciences and Chemical Engineering, Elsevier, 2023.
22. A. Hirsch, Z. Chen and H. Jiao, Angew. Chem., 2000, 39, 3915–3917.
23. P. K. Chattaraj and T. Höltzl, Aromat. Met. Clusters, 2018, 281–306.
24. L. V. Duong, H. T. Pham, N. M. Tam and M. T. Nguyen, Phys. Chem. Chem. Phys., 2014, 16, 19470–19478.
25. H. T. Pham, L. V. Duong and M. T. Nguyen, J. Phys. Chem. C, 2014, 118, 24181–24187.
26. D.-Z. Li, Q. Chen, Y.-B. Wu, H.-G. Lu and S.-D. Li, Phys. Chem. Chem. Phys., 2012, 14, 14769.
27. H.-J. Zhai, Q. Chen, H. Bai, H.-G. Lu, W.-L. Li, S.-D. Li and L.-S. Wang, J. Chem. Phys., 2013, 139, 174301.
28. H. Bai, B. Bai, L. Zhang, W. Huang, H. J. Zhai and S. D. Li, Phys. Chem. Chem. Phys., 2017, 19, 31655–31665.
29. H. Bai, Q. Chen, C. Q. Miao, Y. W. Mu, Y. B. Wu, H. G. Lu, H. J. Zhai and S. D. Li, Phys. Chem. Chem. Phys., 2013, 15, 18872–18880.
30. L. V. Duong and M. T. Nguyen, Phys. Chem. Chem. Phys., 2017, 19, 14913–14918.
31. L. Van Duong, N. N. Tri, N. P. Hung and M. T. Nguyen, J. Phys. Chem. A, 2022, 126, 3101–3109.
32. W. Lu, D. C. H. Do and R. Kinjo, Nat. Commun., 2020, 11, 3370.
33. A. G. Arvanitidis, K. Z. Lim, R. W. A. Havenith and A. Ceulemans, Int. J. Quantum Chem., 2018, 118, 25575.
34. L. Van Duong, D. T. T. Mai, M. P. Pham-Ho and M. T. Nguyen, Phys. Chem. Chem. Phys., 2019, 21, 13030–13039.
35. H.-J. Zhai, B. Kiran, J. Li and L.-S. Wang, Nat. Mater., 2003, 2, 827–833.
36. T. B. Tai, N. M. Tam and M. T. Nguyen, Chem. Phys. Lett., 2012, 530, 71–76.
37. A. P. Sergeeva, D. Y. Zubarev, H.-J. Zhai, A. I. Boldyrev and L.-S. Wang, J. Am. Chem. Soc., 2008, 130, 7244–7246.
38. C.-Y. Gao, Q.-Q. Yan, Q. Chen, Y.-W. Mu and S.-D. Li, J. Cluster Sci., 2024, 35, 1375–1380.
39. M. A. Addicoat and G. F. Metha, J. Comput. Chem., 2009, 30, 57–64.
40. J. H. Holland, Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence, MIT Press, 1992.
41. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
42. Y. A. Jeilani, L. Van Duong, O. M. S. Al Qahtani and M. T. Nguyen, Phys. Chem. Chem. Phys., 2024, 26, 11347–11359.
43. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
44. T. Clark, J. Chandrasekhar, G. W. Spitznagel and P. V. R. Schleyer, J. Comput. Chem., 1983, 4, 294–301.
45. 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, Inc., Wallingford CT.
46. D. Y. Zubarev and A. I. Boldyrev, Phys. Chem. Chem. Phys., 2008, 10, 5207–5217.
47. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
48. L. Rincon, R. Almeida, J. E. Alvarellos, D. Garcia-Aldea, A. Hasmy and C. Gonzalez, J. Chem. Soc., Dalton Trans., 2009, 3328–3333.
49. A. Savin, B. Silvi, F. Coionna and F. Colonna, Can. J. Chem., 1996, 74, 1088–1096.
50. M. Kohout, DGrid, version 5.0, Dresden, 2017.
51. J. C. Santos, W. Tiznado, R. Contreras and P. Fuentealba, J. Chem. Phys., 2004, 120, 1670–1673.
52. L. V. Duong and M. T. Nguyen, Chem. Phys. Lett., 2017, 685, 377–384.
53. G. Monaco, F. F. Summa and R. Zanasi, J. Chem. Inf. Model., 2021, 61, 270–283.
54. S. Jin, B. Chen, X. Kuang, C. Lu, W. Sun, X. Xia and G. L. Gutsev, J. Phys. Chem. C, 2019, 123, 6276–6283.
55. D. T. T. Mai, L. V. Duong, T. B. Tai and M. T. Nguyen, J. Phys. Chem. A, 2016, 120, 3623–3633.
56. D. T. T. Mai, L. V. Duong, M. P. Pham-Ho and M. T. Nguyen, J. Phys. Chem. C, 2020, 124, 6770–6783.
57. R. O. Jones, G. Ganteför, S. Hunsicker and P. Pieperhoff, J. Chem. Phys., 1995, 103, 9549–9562.
58. A. Hirsch, Z. Chen and H. Jiao, Angew. Chem., Int. Ed., 2001, 40, 2834–2838.
59. M. Elian, M. M. L. Chen, D. Michael, P. Mingos and R. Hoffmann, Inorg. Chem., 1976, 15, 1148–1155.
60. K. H. Reddy, D. Usharani, J. F. Nixon and E. D. Jemmis, Chem. – Eur. J., 2011, 17, 9142–9152.
61. C. Heinl, E. Peresypkina, G. Balázs, E. Mädl, A. V. Virovets and M. Scheer, Chem. – Eur. J., 2021, 27, 7542–7548.
62. R. Haunschild and G. Frenking, Mol. Phys., 2009, 107, 911–922.
63. A. S. Ivanov, K. V. Bozhenko and A. I. Boldyrev, J. Chem. Theory Comput., 2012, 8, 135–140.
64. R. G. Parr and R. G. Pearson, J. Am. Chem. Soc., 1983, 105, 7512–7516.
65. R. G. Parr, R. A. Donnelly, M. Levy and W. E. Palke, J. Chem. Phys., 1978, 68, 3801–3807.
66. P. K. Chattaraj, U. Sarkar and D. R. Roy, Chem. Rev., 2006, 106, 2065–2091.
67. A. P. Sergeeva, B. B. Averkiev, H. J. Zhai, A. I. Boldyrev and L. S. Wang, J. Chem. Phys., 2011, 134, 224304.
68. A. P. P. Sergeeva, Z. A. A. Piazza, C. Romanescu, W.-L. L. Li, A. I. I. Boldyrev and L.-S. S. Wang, J. Am. Chem. Soc., 2012, 134, 18065–18073.
69. Z. Chen, C. S. Wannere, C. Corminboeuf, R. Puchta and P. v. R. Schleyer, Chem. Rev., 2005, 105, 3842.
70. C. A. Bauer, A. Hansen and S. Grimme, Chem. – Eur. J., 2017, 23, 6150–6164.
71. R. Nieman, J. R. Carvalho, B. Jayee, A. Hansen, A. J. A. Aquino, M. Kertesz and H. Lischka, Phys. Chem. Chem. Phys., 2023, 25, 27380–27393.
72. A. Pérez-Guardiola, M. E. Sandoval-Salinas, D. Casanova, E. San-Fabián, A. J. Pérez-Jiménez and J. C. Sancho-García, Phys. Chem. Chem. Phys., 2018, 20, 7112–7124.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4dt02251c

---