Searching for double σ- and π-aromaticity in borazine derivatives.

Inspired by the double-aromatic (σ and π) C6H3 +, C6I6 2+, and C6(SePh)6 2+ ring-shaped compounds, herein we theoretically study their borazine derivative analogues. The systems studied are the cation and dications with formulas B3N3H3 +, B3N3Br6 2+, B3N3I6 2+, B3N3(SeH)6 2+, and B3N3(TeH)6 2+. Our DFT calculations indicate that the ring-shaped planar structures of B3N3H3 +, B3N3I6 2+, and B3N3(TeH)6 2+ are more stable in the singlet state, while those of B3N3Br6 2+ and B3N3(SeH)6 2+ prefer the triplet state. Besides, exploration of the potential energy surface shows that the ring-shaped structure is the putative global minimum only for B3N3I6 2+. According to chemical bonding analysis, B3N3H3 +, B3N3I6 2+, and B3N3(TeH)6 2+ have σ and π delocalized bonds. The number of delocalized σ/π electrons is 2/6 for the first, and 10/6 for the second and third, similar to what their carbon analogs exhibit. Finally, the analysis of the magnetically induced current density allows B3N3H3 +, B3N3I6 2+, and B3N3(TeH)6 2+ to be classified as strongly σ aromatic, and poorly π aromatic compounds. This journal is © The Royal Society of Chemistry.

Introduction

In 1979, Schleyer and co-workers introduced the concept of double aromaticity. They showed that the 3,5-dehydrophenyl cation (1, Scheme 1) exhibits two orthogonal Hückel frameworks, consisting of 6π and 2σ electrons.[1] Since then, many compounds exhibiting double aromaticity have been predicted theoretically,[2-11] with some even being experimentally characterized.[12-17] For instance, in 1988, Martin and Sagl[15] reported the synthesis of the hexaiodobenzene dication, C6I62+ (2, Scheme 1), a double-aromatic compound with two concentric delocalization circuits, consisting of 6π electrons delocalized around the C6 ring, and 10σ electrons delocalized around the I6 ring. This finding paved the way to synthesize and characterize, by X-ray diffraction, the first double-aromatic “bench-stable” molecule, the hexakis(phenylselenyl)benzene dication[16] (3, Scheme 1).

Scheme 1

Structures of double-aromatic benzene derivatives (1–3) and their borazine derivative analogues. Dashed lines represent the σ delocalization circuits.

Here we evaluate the structure, stability, chemical bonding, and magnetic response properties of B3N3H3+, B3N3Br62+, B3N3I62+, B3N3(SeH)62+ and B3N3(TeH)62+ (4, 5, 6, 7 and 8, Scheme 1) systems in pursuit of expanding the family of double-aromatic compounds.

As we can see in Scheme 1, compounds 4–8 are analogues to the double-aromatic benzene derivatives referred above (1, 2 and 3, Scheme 1), where C6 is replaced by the B3N3 ring. It is important to note that the aromatic character of borazine has been a much-debated topic in the past with contrasting conclusions.[18-27] There is even disagreement among studies based on the same criteria, i.e. magnetic. For instance, various studies have evaluated current density maps, both global and dissected (σ and π), and concluded that borazine does not support a net diatropic current (aromaticity indicator) around the B3N3 ring. Instead, it displays three local π delocalization pathways (around the N atoms).[18-20] Given the large difference in electronegativity between the two participating elements, Δχ(N,C) = 1.0,[28] it sounds reasonable that the π electrons are more polarized to the N atoms. However, analysis of both the induced magnetic field[26,27] and the intensity of the ring current strengths (RCS)[23,24,26] suggests that borazine does sustain a net diatropic ring current (less intense than that of benzene). Given these discrepancies, why would we assume that the compounds 4–8, in this work, would have a double aromatic character? Our answer is based on the quantitative analyses on the aromaticity of borazine available in the literature. The most modern analyses, based on the energetic criteria (aromatic stabilization energy, ASE), agree that the aromaticity of borazine is at least 20% benzene's.[21,22,26] For instance, Baranac-Stojanović et al. reported extra cyclic resonance energy (ECRE) values, at B3LYP/6-311+G(d,p) level, of 89.35 and 20.18 kcal mol−1 for benzene and borazine, respectively.[22] Furthermore, in 2013, our group estimated that the aromaticity of borazine is about 30% benzene's, according to the energetic and magnetic criteria.[26] The first one, using isodesmic reactions to estimate ASEs, and the second one, through the analysis of three properties: the induced magnetic field, the nucleus independent chemical shift (NICS) and the RCS. These results contrast slightly with RCSs of benzene and borazine (12.8 and 2.1 nA T−1), calculated at RHF level, by Monaco et al.[29] However, they agree with both RCS values for benzene and borazine reported most recently by Monaco et al. (11.8 and 3.0 nA T−1)[23] and with our current calculations (11.8 and 3.6 nA T−1). Under this energetic and magnetic evidences, we consider that the B3N3 ring, of the systems studied here, should sustain at least a slight π-aromatic character. The latter, together with the fact that the chemistry of borazine, and its derivatives, is well established,[30-32] supports our proposal of compounds 4–8, as potential aromatic derivatives of borazine.

Computational details

In order to evaluate the stability, chemical bonding and aromatic character of the proposed systems, we used the following methodologies:

(a) The potential energy surface (PES) of systems, 4, 6 and 8, was explored via the AUTOMATON program.[33] This program uses a probabilistic automaton method to generate an initial population that evolves to the best individual (global minimum, GM) through genetic operations. The searches were performed via the PBE0 [34] functional in conjunction with the SDDALL[35,36] basis. The low-lying structures were subsequently re-minimized at the PBE0-D3 [37]/def2-TZVP[38] level. We used an initial population of 5X (X = number of atoms of the system), this small population showed to be enough to successfully identify the lowest energy structures of many clusters and molecules.[39-43] Additionally, we evaluated the dynamic behavior of compounds 4, 6 and 8 (for an insight into their kinetic stability), employing Born–Oppenheimer molecular dynamic (BOMD) simulations.[44] Geometry optimizations, vibrational frequency and BOMD calculations were carried out using the Gaussian 16 program.[45]

(b) Chemical bonding was analyzed through the adaptive natural density partitioning (AdNDP) method,[46,47] which is an extension of the natural bond orbital (NBO) analysis.[48-50] AdNDP allows for electronic density partitioning of a molecular system in terms of n-center two-electron (nc-2e) bonds, with n ranging from one to the total number of atoms in the molecule. Thus, AdNDP recovers the Lewis elements (lone pairs, and 2c-2e bonds), as well as the delocalized nc-2e bonds. Therefore, the number of nc-2e (n > 2) delocalized bonds allows to classify a system as aromatic (4N + 2) or anti-aromatic (4N), following Hückel's rule.[51-53]

(c) The aromatic character was also assessed through the magnetic criteria. Current density maps, induced by an external magnetic field, were calculated for this purpose. These maps allowed us to easily observe the diatropicity/paratropicity of the examined systems, and to correlate these magnetic patterns with their aromatic/antiaromatic character.[54-56] Additionally, the ring current flow was integrated to estimate the ring current strength (RCS) and to use it as a quantitative criterion for aromaticity.[57-59] In general, a positive and a negative RCS value indicate an aromatic and an anti-aromatic character, respectively. Similarly, the larger the absolute value of the RCS, the higher the (anti) aromatic character of the studied molecule. However, it is more difficult to establish an exact RCS value for nonaromatic systems. The value must be close to zero, but how close? To answer this, we need to compare them with suitable references. In this paper, we have calculated the RCSs of a series of five-member C4XH4 rings, for which Mo and Schleyer reported their “extra cyclic resonance energy” (ECRE) to characterize and measure the extra stabilization (aromaticity).[60] ECRE was computed using block-localized wave function (BLW) method[61] and adequate isodesmic reactions. The simple linear correlation of the RCSs vs. ECREs is excellent (r2 = 0.97). Furthermore, using the linear equation of this correlation, a value of 0 kcal mol−1 ECRS (nonaromatic) would correspond to an RCS value of 0.88 nA T−1 (see Table S1 and Scheme S1 in the ESI†). This analysis allows us to consider RCS values higher than 1.5 nA T−1 here, as indicative of aromatic character. In addition, our group has shown that RCS correlates very well with NICS, when the latter are measured properly.[62-65]

Results and discussion

We began by analyzing compound 4, which is not the global minimum on the B3N3H3+ PES, but it lies at 49.3 kcal mol−1 above the lowest singlet state energy structure at PBE0 [34]-D3 [37]/def2TZVP level.[38] However, 4 is a true minimum on its PES since all its vibrational frequencies are positive. The 22 lowest singlet state energy structures of B3N3H3+, identified by the AUTOMATON program, are reported in Fig. S1.† It is important to note that AUTOMATON did not identify structure 4, since this program focuses on identifying the GM structure. However, AUTOMATON allowed discarding 4 as the putative GM.

Interestingly, the structure of 4 is similar to that of its analog 1. As with the carbon angles bearing no hydrogens at 1, the angle widens in the boron of compound 4. The optimized structure at PBE0-D3/def2TZVP level shows angles NBN to be 167°. The B–B Wiberg bond index (WBI)[66,67] (0.45|e| at PBE0-D3/def2TZVP level) is similar, to that obtained for the dehydrogenated carbon atoms in 1, which according to Schleyer, has a σ 3c-2e delocalized bond. This result supports a possible in-plane 3c-2e σ-bonding interaction involving the B3 ring.

According to the AdNDP results reported in Fig. 1, the 26 valence electrons (13 pairs) are distributed in 6 classical B–N σ-bonds (2c-2e) and 3 classical N–H σ-bonds. All other bonds are delocalized. There are three 6c-2e π-bonds (involving the B3N3 ring), and one 3c-2e σ-bond (involving the B3 ring). These results agree with the WBI predictions. Note that 6π/2σ electron counting all conform to the (4n + 2) Hückel rule, providing further support for the double-aromatic character of 4. Interestingly, AdNDP predicts that chemical bonding in compound 4 is similar to that in 1, the 3,5 dehydrophenyl cation (C6H3+; Fig. 1b).

Fig. 1

AdNDP bonding pattern of (a) B3N3H3+ (4) as compared to that of (b) C6H3+ (1). Occupation numbers (ONs) are shown.

Fig. 2 depicts the magnetically induced current densities calculated in planes located at 1 a0 (1 bohr) above the molecular plane of 4 (1 is also included for comparison). The ring current strength (RCS) was computed by integrating the ring current flux that passes along interatomic surfaces defined by the quantum theory of atoms in molecules (QTAIM).[58,59] Diatropic currents are assumed to circle clockwise, and the paratropic ones to circle anticlockwise. The current density map of 4 exhibits two diatropic ring currents (one σ and one π) distributed around the molecular ring (Fig. 2a and b). As expected, compound 1 also exhibits diatropic σ and π ring currents (Fig. 2c and d) with a σ-current flow distribution pattern similar to 4. However, π ring current flow patterns are different for 1 and 4: the ring current is continuous and clearly noticeable at 1, but at 4, there are two types of current, the diatropic local one on each nitrogen atom, and the ring current. Since the local is stronger than the ring current, the pattern for 4 does not exhibit the same continuity as in 1. How intense are these current flows, and therefore, how aromatic is system 4? The σ-RCS/π-RCS value of 5.4/1.5 nA T−1 predicts system 4 to be a σ-aromatic/poorly π-aromatic compound; whereas σ-RCS/π-RCS value of 11.1/9.9 nA T−1 predicts system 1 to be a σ-aromatic/π-aromatic compound. This assignment is made on the basis of the reference σ-RCS/π-RCS values of 0.5/11.7 nA T−1 for benzene, which is a σ-nonaromatic/strongly π-aromatic system. It is important to note that 1 is the global minimum structure for the C6H3+ combination (at the CCSD(T)-F12b/cc-pVTZ + ZPE(B3LYP/cc-pVTZ)//CCSD(T)/cc-pVTZ level).[68] Thus, in the light of our current density analysis, it is possible to suggest that the high σ and π aromatic character must be contributing to ring shaped stabilization of C6H3+. Whereas system 4, despite its substantial σ-diatropic character, is not an energetically competitive isomer.

Fig. 2

Vector maps of dissected magnetically induced current density (MICD) at 1 a0 above the molecular plane of (4) B3N3H3+ (top) and (1) C6H3+ (bottom).

We now focus on discussing our results for the set of persubstituted borazine dications, B3N3Br62+ (5), B3N3I62+ (6), B3N3(SeH)62+ (7) and B3N3(TeH)62+ (8). From this series, only structures 6 and 8 are more stable in the singlet state, while structures 5 and 7 prefer the triplet state (at PBE0-D3/def2TZVP level). Therefore, 5 and 7 are discarded as suitable candidates for double-aromatic systems.

After exploring its PES, AUTOMATON predicts 6 as the lowest energy structure for the B3N3I62+ combination, whereas 8 lies 68.7 kcal mol−1 above the lowest energy singlet-state structure for B3N3(TeH)62+ combination (all structures identified by AUTOMATON within 70 kcal mol−1 above the global minimum for B3N3(TeH)62+ are reported in Fig. S2†). It is important to note that canonical AUTOMATON procedure provided dissociated structures as the lowest energy isomers for these systems. Thus, we performed a guided search, which consists in replacing the H by I or by the Te–H fragment in the 60 lowest energy structures identified by AUTOMATON for borazine, subsequently the resulting structures were used as initial population in the AUTOMATON's PES exploration of B3N3I62+ and B3N3(TeH)62+. The second isomer of B3N3I62+ is the triplet-state ring shaped C2v structure and lies 2.0 kcal mol−1 above 6. Singlet and triplet isomers under 30 kcal mol−1 identified by AUTOMATON for B3N3I62+ combination are reported in Fig. S3 and S4 in the ESI,† respectively.

At the optimized structures (at PBE0-D3/def2TZVP level) of 6 and 8, the B–N distances are 1.42 and 1.43 Å respectively, whereas I–I, and Te–Te distances are 3.49 and 3.55 Å, respectively. The WBI value for B–N in compound 6 and 8 are 0.94, very similar to the B–N WBI value in borazine, 1.02, suggesting that B3N3 ring in compound 6 and 8 is bonded in a similar way to the B3N3 ring of the borazine. On the other hand, WBI values for I–I and Te–Te bonds are both 0.13, suggesting a covalent bonding interaction between these atoms.

According to AdNDP, the bonding pattern of 6 is composed of localized and delocalized elements (Fig. 3a). The localized part is found as twelve lone pairs in the I atoms (two LPs on each one), six classical B–N 2c-2e σ-bonds, three 2c-2e B–I σ-bonds, and three 2c-2e N–I σ-bonds. The delocalized elements are found as three 6c-2e π-bonds (around the B3N3 ring), and five 6c-2e σ-bonds (around the I6 ring). Thus, AdNDP analysis classifies 6 as a 6π/10σ electron system, doubly satisfying the (4n + 2) Hückel's rule, providing further support for its double-aromatic character. Note that AdNDP predicts that chemical bonding in compound 6 is similar to that in 2 (C6I62+; Fig. 3b). Moreover, under the same analysis scheme, system 8 shows a very similar chemical bonding pattern to the one shown for system 6, as can be seen in Fig. S5 in the ESI.†

Fig. 3

AdNDP bonding pattern of (a) B3N3I62+ (6) as compared to that of (b) C6I62+ (2). Occupation numbers (ONs) are shown.

Fig. 4 depicts the magnetically induced current densities calculated in planes located at 1 a0 (1 bohr) above the molecular plane of 6 (2 is also included for comparison). The ring current strength (RCS) was computed by integrating the ring current flux that passes along C–C, B–N and I–I interatomic surfaces, defined by the quantum theory of atoms in molecules (QTAIM). The current density maps of 6 exhibit two diatropic ring currents, one outer-σ (o-σ) and intense distributed around I6 ring, and the other inner-π (i-π) and weak distributed around B3N3 ring (Fig. 4a and b). Compound 2 exhibits a similar o-σ-ring current flow pattern to 6 but a different i-π ring current flow pattern, which is strong in 2. Moreover, the o-σ-RCS/i-π-RCS value of 14.3/3.4 nA T−1 predicts system 6 to be a σ-aromatic/poorly π-aromatic compound, whereas o-σ-RCS/i-π-RCS value of 14.9/10.5 nA T−1 predicts system 2 as a σ-aromatic/π-aromatic compound. For compound 8, a similar pattern to 6 is found, with o-σ-RCS/i-π-RCS value of 8.5/3.6 nA T−1 (see Fig. S6 in the ESI†). Thus, according to ring currents, the borazine derivatives 6 and 8 are σ-aromatic and weakly π-aromatic systems.

Fig. 4

Vector maps of dissected magnetically induced current density (MICD) at 1 a0 above the molecular plane of B3N3I62+ (top) and C6I62+ (bottom).

Finally, BOMD simulations were performed to evaluate the dynamic behavior of compounds 4, 6 and 8 (Movies S1, S2 and S3 in ESI†). The corresponding movies are available in the ESI.† We found that a cyclic arrangement persists over the period of time considered for the dynamics (30 ps), at 900 K. Interestingly, the GM for the combination B3N3(TeH)62+ is not kinetically stable, dissociating in the first steps of the BOMD at 900 K (see Movie S4 in the ESI†). These results suggest good kinetic stability for these compounds, which in conjunction with the evidence that 6 is the putative global minimum for the B3N3I62+ combination, hold promise for its experimental realization.

Conclusions

In summary, we examined five double σ and π aromatic candidates with ring-shaped structures, B3N3H3+ (4), B3N3Br62+(5), B3N3I62+(6), B3N3(SeH)62+ (7) and B3N3(TeH)62+ (8), derived by isoelectronic substitution of the C6 ring of the parent species, C6H3+, C6Br62+, C6I62+ C6(SeH)62+, and C6(TeH)62+. From this series, only 4, 6, and 8 structures are more stable in singlet state, whereas 5 and 7 prefer the triplet state. However, the energy of 4 and 8 is 63.2 kcal mol−1 and 68.7 kcal mol−1 higher than that of the corresponding lowest energy singlet state structure respectively, but structure 6 is the putative global minima on its PES. Moreover, the ab initio BOMD simulations revealed that 6 and 8 structures are rather rigid, persisting over the period of time considered for the dynamics (30 ps), at 900 K. Hence, these results suggest that it might be possible to detect 6 as the thermodynamically and kinetically viable structures in the gas phase. Another promising approach would be to replace TeH by TeR ligands, to include steric protection by using bulky R-groups, considering that 8 shows reasonable kinetic stability (according to the BOMD simulations).

Conflicts of interest

There are no conflicts to declare.