Dariusz W.
Szczepanik
*a,
Marcin
Andrzejak
a,
Justyna
Dominikowska
b,
Barbara
Pawełek
c,
Tadeusz M.
Krygowski
d,
Halina
Szatylowicz
e and
Miquel
Solà
f
aK. Gumiński Department of Theoretical Chemistry, Faculty of Chemistry, Jagiellonian University, Gronostajowa 2, 30-387 Kraków, Poland. E-mail: dariusz.szczepanik@uj.edu.pl
bDepartment of Theoretical and Structural Chemistry, Faculty of Chemistry, University of Lodz, Pomorska 163/165, 90-236 Łódź, Poland
cDepartment of Plant Cytology and Embryology, Institute of Botany, Jagiellonian University, Gronostajowa 9, 30-387 Cracow, Poland
dDepartment of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warszawa, Poland
eFaculty of Chemistry, Warsaw University of Technology, Noakowskiego 3, 00-664 Warszawa, Poland
fInstitut de Quìmica Computacional i Catàlisi, Universitat de Girona, C/Maria Aurelia Capmany, 69, 17003 Girona, Catalonia, Spain
First published on 17th October 2017
In this study the recently developed electron density of delocalized bonds (EDDB) is used to define a new measure of aromaticity in molecular rings. The relationships between bond-length alternation, electron delocalization and diatropicity of the induced ring current are investigated for a test set of representative molecular rings by means of correlation and principal component analyses involving the most popular aromaticity descriptors based on structural, electronic, and magnetic criteria. Additionally, a qualitative comparison is made between EDDB and the magnetically induced ring-current density maps from the ipsocentric approach for a series of linear acenes. Special emphasis is given to the comparative study of the description of cyclic delocalization of electrons in a wide range of organic aromatics in terms of the kekulean multicenter index KMCI and the newly proposed EDDBk index.
Among all the ways to quantify aromaticity through the “ground-state” criteria (energetic, structural, and electronic), there are descriptors that have gained enormous recognition and wide acceptance: the aromatic stabilization energy (ASE),2 harmonic oscillator model of aromaticity (HOMA),11–13 and different types of the multicenter delocalization indices (MCDI).14–18 In turn, one of the most popular aromaticity descriptors based on the response properties is the nucleus-independent chemical shift (NICS) and its various derivatives.4,19,20 Despite the unquestionable success and popularity of these aromaticity indices, some of their imperfections still cuts back their applicability to relatively small and simple systems. In particular, design of isodesmic and homodesmotic reaction scenarios to determine ASE is very difficult in practice and it leaves room to a lot of arbitrariness.2 The principal problem with HOMA, in turn, is the necessity of parametrization of bond lengths for an idealized reference molecule, which obviously cannot be chosen unambiguously. Consequently, the practical use of HOMA is limited to aromatic and heteroaromatic systems since the parameters for chemical bonds with metal atoms are not available.3 Furthermore, the parametrization of HOMA should be performed using the same quantum-chemical method as used in calculations of equilibrium geometries of the molecule under study, since routinely computed HOMA with the experimentally determined parameters is bound to suffer from large and unsystematic errors.13 The magnetic-based measures of aromaticity, like NICS, have also been criticised for their complexity (NICS relies on the condensation of potentially complicated patterns of induced currents to a single number, especially in the case of fused aromatic rings),21–26 methodological shortcomings (aromaticity evaluation using NICS is limited mainly to planar units of similar size),6,20,27 and interpretative mistiness (quoting Prof. P. Bultinck, “…a multicenter delocalization is a necessary condition for a diatropic ring current to exist, provided that there are proper virtual molecular orbitals to excite to.”).28–32 Finally, the DI-based indices (especially MCDI), unlike the aromaticity measures described above, enable one to study most of the types of aromaticity that can be found in literature.6 Also, the MCDIs are the only descriptors that passed a set of rigorous tests for aromaticity quantifiers.33–35 Unfortunately, the main disadvantage connected with calculation of MCDI is their computational cost – at the level of one-determinant wavefunction with moderate basis set one cannot use multicenter delocalization indices to evaluate aromaticity of the molecular fragments containing more than a dozen of atoms.6 Thus, for large molecular rings, the analysis of aromatic stabilization using MCDI is applicable only within the framework of simple approximations like Hückel's or pseudo-π methods and for limited range of cases.36
Recently the original method of the electron density of delocalized bonds (EDDB) has been proposed to facilitate quick qualitative analysis of different bond-conjugation patterns and to provide a bird's-eye view on the global aromaticity and resonance effects in molecular systems that due to their size and complex structure are the major challenge for the currently used tools.37,38 However, in a number of preliminary tests EDDB turned out to be highly capable of providing also a quantitative evaluation of electron delocalization in many diversified aromatic rings and simultaneously free from the aforementioned shortcomings of the commonly used aromaticity indices.39 This work provides a comprehensive view on the performance of the electron density of delocalized bonds in quantification of local aromaticity by comparative study with other well-known descriptors based on structural, magnetic and electronic criteria. Special emphasis is given to the analysis of cyclic delocalization of π-electrons in terms of the newly proposed EDDB-based index and the kekulean multicenter index. Since the former does not strictly take into account the cooperativity of all atomic centers in cyclic delocalization of electrons, as the latter does, the results will show if and to what extent the multicenter sharing effects are important for reliable description of organic aromatics.
ED(r) = EDLA(r) + EDLB(r) + EDDB(r). | (1) |
(2) |
(3) |
Ω = {(A1 − A2),(A2 − A3),(A3 − A4),(A4 − A5),(A5 − A1)}. | (4) |
Fig. 1 (a) Different levels of electron delocalization from the one-electron density decomposition scheme by eqn (1) for pyridine (only valence electrons); (b) global EDDB(r) function for the LEU-LYS-GLU-GLN-PRO-ARG-HIS-PHE-TYR-TRP decapeptide (colored fragments indicate aromatic rings). Method: HF/6-311G**//PM6. |
It has to be emphasized that EDDB is far more efficient than MCDI, especially in the case of highly accurate wavefunctions of large molecular systems. For instance, the calculation of MCDI takes from a dozen of seconds to several hours depending on the size of the ring and the computational method used, while the EDDB calculation takes less than 1 s for all the aromatic rings considered in this study regardless of their size.39 Even the determination of global aromaticity/resonance effects in molecules containing hundreds of atoms, like the decapeptide depicted in Fig. 1b, is very fast and takes less than 40 s if the appropriate threshold for the 2cBO occupations is used and much less than 10 minutes otherwise.39 Such speed-up is possible because within the EDDB formalism the multicenter electron sharing is approximated by means of decoupled three-atomic local resonances representing conjugations between the adjacent bonds only (cf. the Bridgeman–Empson method).60 Further reduction of computational time can be achieved by using the natural minimal basis (NMB) instead of the default (full) NAO-representation.
Fig. 2 The T1 test set of molecular rings used to benchmark the performance of different aromaticity descriptors. |
We have chosen the following descriptors as the most representative and commonly used aromaticity measures that can be found in literature:
• Harmonic oscillator model of aromaticity (HOMA),11–13 which is a normalized measure of the energetic consequences of deviations of bond lengths in the molecular ring from the corresponding optimum values for an idealized aromatic system.
• The axial component of the nucleus-independent chemical shift calculated at 1 Å above the ring centroid, NICS(1)zz,4,19,20 which quantifies diatropicity/paratropicity of the induced ring current by means of the effective magnetic shielding.
• The kekulean multicenter index (KMCI),14–18 which directly quantifies cyclic delocalization of electrons in aromatic rings. Our preliminary studies indicate that the KMCI values for all aromatic species from both test sets tightly correlate (R2 = 1.00) with those of MCI – the multicenter index originally proposed by Bultinck et al.15 that implicitly takes into account the cross-ring delocalization of electrons. However, since MCI is far more computationally expensive and offers virtually no advantage over KMCI (at least for the studied systems) its use in this benchmark is not necessary; the MCI values are included in Table S1 in ESI.† In view of the well-known problems of the multicenter indices with the ring-size extensivity, in order to perform the correlation and principal component analysis (PCA) involving the entire T1 set of 5-, 6-, and 7-MR molecules we used the nth root of the original index, denoted by KMCI1/n (see Fig. 3a and 4a), as suggested in the literature.6,16
• Average two-center index (ATI),15 which measures the average effect of the electron delocalization between three para-related atoms in a benzenoid unit. Its use is thus limited to 6-membered rings only. ATI has the same theoretical foundations as the well-known para-delocalization index (PDI) by Poater et al.5,62 but it involves a Hilbert-space partitioning within the basis of atomic orbitals instead of Bader's atoms-in-molecule (AIM)65 or the fuzzy-atomic space (FAS).61
• Shannon aromaticity (SA),63 which measures the Kullback–Leibler distance of the bonding electron density distribution in aromatic ring from uniformity; SA is a non-referential index.
• Fluctuation index of aromaticity (FLU),5,64 which has a similar interpretation as the SA index but depends upon parameters determined for an idealized aromatic system.64
• Ellipticity index of aromaticity (EL),66 which measures bonding electron density deformations reflected in negative eigenvalues of the Hessian matrix of the electronic density in the bond critical point (BCP).
• Density of the total electron energy at the ring critical point (HRCP),67 which is one of the AIM parameters that has been proven to serve as a quantitative measure of aromaticity in a number of molecular rings.68
Formal definitions of all the above listed indices can be found in the source papers as well as in the manual of the MultiWFN program,69 which has been used to calculate HOMA (with parameters for the C–C bond calculated consistently according to ref. 9) as well as ATI and KMCI (both in the NAO basis). To calculate HRCP, FLU, SA, and EL the analysis of electron density distribution within the framework of quantum theory of atoms in molecules has been performed using AIMAll program.70 Although the electronic indices calculated using different partitioning schemes may sometimes give rise to different aromaticity predictions,71 our preliminary studies show that the NAO-based indices like KMCI and ATI are fully equivalent to their AIM- and FAS-based counterparts, i.e. IRing and PDI, respectively (the r-squared coefficients are very close to 1.00 in all cases); the IRing and PDI values are included in Table S1 in ESI.† The NICS(1)zz values have been determined for all systems (using the Gaussian 09 package72) except for polycyclic aromatic hydrocarbons, for which NICS is ill-defined, owing to its proven non-local character caused by mixing of different ring currents.21–25 The EDDBk indices have been obtained using the NBO6 software73 and the script program written by one of the authors (DS).74 The corresponding EDDB(r) maps have been generated by means of the standard tools from the Gaussian 09 package (formchk and cubegen programs). The CAM-B3LYP75/def2-TZVPP76 calculations with full geometry optimizations have been performed using Gaussian 09. The correlation and principal component analyses based on appropriately scaled aromaticity indices (scaling factors can be found in Table S1 in ESI†) have been carried out using the R-package77 and their results are summarized in Fig. 3 and 4, while the complete results of the benchmark calculations are collected in Tables S1 and S2 in ESI.†
Test set T2 contains more diverse types of aromatic systems and has been used to assess if and to what extent the cyclic delocalization patterns in aromatic rings predicted by KMCI can be reproduced by the newly proposed index EDDBk; a more comprehensive study involving the AIM-based counterpart of KMCI and other indices based on different criteria of aromaticity can be found elsewhere.33–35 The T2 set is based on the collection of tests originally designed by Solà et al.33 to evaluate new aromaticity indicators proposed in the literature and it accumulates chemical experience about the expected trends in aromaticity changes in the following systems: distorted benzene, substituted benzene, metal complexes, penta- and heptafulvenes, claromatic systems, heteroaromatic systems, as well as the aromatic transition states in selected chemical reactions. All the test molecules are depicted in Fig. 6–8. KMCI and EDDBk indices were calculated at the B3LYP/6-311++G** theory level (equilibrium geometries) using Gaussian 09 and are listed in Fig. 6 and 7 (results for distorted benzene and two chemical reactions); Fig. 8 presents a summary of the results for other species from the test set, which are displayed and briefly described in Fig. S1–S5 in ESI.†
Even a cursory glance at Fig. 3a (as well as Fig. 4b and g) indicates that the ring-size extensivity issue6 makes HRCP and NICS(1)zz incomparable with the rest of aromaticity descriptors when molecular rings of different size are considered.78 On the other hand, reasonable linear correlation (with R2 coefficients close to or greater than 0.90) can be observed within the set of SA, EL, FLU and HOMA indices, while their correlation with EDDBk is mainly non-linear in character, but still relatively tight (up to R2 = 0.96 for HOMA). This is clearly shown in Fig. 4c–f. The fact that aromaticity changes predicted by different aromaticity measures may not necessarily be linearly proportional to each other is quite obvious. It should be noticed, however, that here the exponential relation between EDDBk (a quantity based on bond-order orbitals) and HOMA (an index involving bond lengths) explicitly refers to the bond-distance/bond-order correlation originally established by Pauling.79 As follows from Fig. 3a, in turn, the renormalized multicenter index, KMCI1/n, seems to fall slightly behind EDDBk, SA, EL, FLU, and HOMA with r-squared coefficients in the range of 0.70 to 0.90 (excluding HRCP and NICS(1)zz), but it still performs dramatically better than the original KMCI for size-differentiated aromatic rings (see Fig. 4a and Table S2a in ESI†).
The results of PCA indicate that PVE1C for the entire set of aromaticity indices calculated for all aromatic rings from T1 equals 78% (see Fig. 3a). However, consecutive eliminations of indices with the smallest contributions to the first component (and significant contributions to the second one) from the set systematically improve PVE1C. In particular, by eliminating the HRCP measure from the set PVE1C grows to a value of 87%, then by exclusion of NICS(1)zz it reaches 91%, next, elimination of EL increases PVE1C up to 93%, and so on; red numbers indicate PVE1Cs for the last two indices remaining after exclusion of all other aromaticity measures (here FLU and SA). In fact, it should not be surprising that two mutually linked indices (in a sense both quantify uniformity of the electron distribution in aromatic ring) remain at the end of the procedure. But for our purposes the final effect is not as important as the partial result of elimination of a specific descriptor: the smaller the PVE1C growth, the better. In this context, EDDBk (eliminated as the fourth in order) gives rise to relatively small increase of PVE1C by 2 percent points (pp), i.e. from 93% to 95%; basically the same is true for the T1 subsets of 5- (+2pp), 6- (+2pp) and 7-membered rings (0pp). Generally, if one regards size-differentiated rings separately (Fig. 3b–d) PVE1Cs for the entire set of aromaticity indices are close to 90%, which means that all the criteria of aromaticity used in this study give a more or less consistent picture of aromatic stabilization for the set T1. Moreover, a closer inspection of the antecedent PVE1C values reveals that in each case EDDBk is in the subset of only a few descriptors with the rate of at least 96% of variance explained in the case of monocyclic aromatics (Fig. 3b and d) and at least 95% in the polycyclic aromatic hydrocarbons (Fig. 3c).
The effect of the ring-size extensivity issue of the original (unnormalized) kekulean multicenter index is clearly shown in Fig. 4a, in which the correlation with the EDDBk index is analyzed. Apparently, KMCI1/n is superior to KMCI (R2-values in brackets) if one wants to compare aromaticity of rings of different size; at the same time renormalization of KMCI seems to give only a slight adjustment to the r-squared coefficients when we consider rings of different size separately. One should realize, however, that the situation dramatically changes e.g. in the case of NICS(1)zz, where the correlation coefficients for the entire test set T1 speaks in favor of the original KMCI. Indeed, accordingly to Table S2a in ESI,† renormalization of KMCI decreases R2 from 0.67 (exponential correlation) to 0.36 (linear correlation). This particular should not be surprising, since the nucleus-independent chemical shift is known to share the lack of ring-size extensivity with the unnormalized multicenter index.6,20 In contrast, when regarding 5- and 7-MRs separately renormalization of KMCI has negligible impact on the correlation with NICS(1)zz (see Tables S2b and d in ESI†). It should also be mentioned that excluding 1,3-cyclopentadiene from the subset of 5-MRs significantly improves the correlation with NICS(1)zz of both cyclic delocalization measures, KMCI and EDDBk, i.e. the corresponding correlation coefficients rise from 0.79 to 0.97 (KMCI) and from 0.86 to 0.94 (EDDBk), i.e. by 23% and 9%, respectively. When the AIM-based counterpart of KMCI (i.e. the IRing index) is considered, the R2 value increases by nearly 60% (unpublished results). Thus, it seems that EDDBk performs slightly better than KMCI regarding the correspondence between magnetic and electronic criteria of aromaticity, especially when both Hückel's aromatic and antiaromatic systems are taken into account.
A closer look at Fig. 4e and f reveals that the correlation of EDDBk with FLU and HOMA is relatively tight regardless of whether we consider rings of different size separately or collectively. Additionally, in the case of SA (Fig. 4c), excluding the central ring (27) in triphenylene from the 6-MR subset of T1 significantly improves the correlation with EDDBk (as well as with other indices like HOMA and FLU) giving rise to R2 = 0.95. On the other hand, it follows from Fig. 4a, b, d and h that in the case of polycyclic aromatic hydrocarbons the correlation with KMCI1/n, HRCP, EL and ATI is significantly weaker. This effect is associated with a more general problem of the definition of local aromaticity in polycyclic systems.21,22 A good example here are linear acenes as they have been a subject of heated debate in literature, owing to dramatic discrepancies between local aromaticity descriptions provided by different aromaticity criteria (“the anthracene problem”).21,22,80–83
The aromaticity indices based on magnetic and structural criteria predict increasing aromaticity going from terminal to the central ring, while the electronic aromaticity indices like IRing, PDI, or HRCP predict the opposite trend.33,34,84 It has recently been demonstrated, however, that in the latter group of indices the results dramatically depend on the choice of the exchange–correlation functional at the DFT theory level and as such they are somewhat less reliable at least in the case of polyacenes.84 In this context, the EDDBk index gives virtually the same predictions as HOMA and FLU (regardless of the level of the theory). Furthermore, even the qualitative comparison of the global EDDB(r) maps with the magnetically induced ring-current density maps from the ipsocentric approach4,23,85 presented in Fig. 5 clearly show that, despite fundamental methodological differences (induced ring-current density is a response property that occurs only in the presence of an external magnetic field while the cyclic delocalization of electrons is a ground-state property that is present irrespective of the presence or absence of an external magnetic field), both approaches lead to the same conclusions about local aromaticity of acenes as those from HOMA and FLU calculations.
Fig. 5 Contours of the (global) electron density of delocalized bonds and the corresponding maps of the ring-current density from the ipsocentric approach for the first four acenes (only the π-type molecular orbitals were taken into account). Method: B3LYP/6-311G**, equilibrium geometries. CRD(r) maps reproduced from ref. 85 with permission from Elsevier. |
Fig. 6 Different benzene distortions together with the corresponding KMCI and EDDBk values as well as the EDDB(r) contour maps. KMCI values have been multiplied by 103. Method: B3LYP/6-311++G**. |
Fig. 7 presents plots of KMCI and EDDBkvs. the reaction coordinate for both reactions. Analysis of these plots leads to the conclusion that despite some minor differences in the shape of curves, both aromaticity indices precisely identify the σ-aromatic transition state and the π-aromatic product of the acetylene trimerization (here the product is correctly identify as more aromatically stabilized than the corresponding TS) as well as the aromatic TS in the Diels–Alder cycloaddition. The EDDB(r) contour maps show very clearly that TS in the DA reaction has a boat-like structure with characteristic cyclic delocalization pattern in the plane between butadiene and ethylene fragments.88 Moreover, although both TS and the product of the trimerization reaction are to a similar extent stabilized by resonance (according to the EDDBk values, the difference is less than 10%) it is quite obvious even from the first look at the EDDB(r) contours that TS is stabilized by cyclic delocalization of σ-electrons while the final product represents a typical for organic aromatics π-delocalization pattern.
Three additional groups of aromatic systems complete the T2 set of molecules: heteroaromatics, claroaromatics and fulvenes. The first one was originally proposed to predict the proper trend of aromaticity changes along a well-established heteroaromatic series including five aza-derivatives of benzene and five heterocyclic compounds of type C4H4X (where X = CH−, NH, O, CH2, BH, CH+); the second group was designed to test the effect of fusing aromatic rings represented by five Clar systems, while the third subset was used to assess the expected trend of aromaticity in 5-MR and 7-MR fulvenes with different substituents.33 In general, all the electronic aromaticity indices were reported to pass the three tests.33,34 The only exceptions concern the series of pentafulvenes, which are particularly difficult to assess by aromaticity quantifiers as they display a tunable aromatic character (sometimes being termed “aromatic chameleons”92). Magnetic and structural indices are in most cases also in line with the expected trends for all three tests, except a single incorrect prediction that C4H4NH is more aromatic than C5H5−.33,34
How does the EDDBk index deal with the above mentioned tests? The answer to this question is provided by Fig. 8, which displays a summary of calculations for all 47 test systems used to compare the performance of EDDBk and KMCI; detailed results with comments are collected in Fig. S1–S5 in ESI.† The presented results leave no doubt: the index based on the electron density of delocalized bonds predicts exactly the same trends of aromaticity changes as the kekulean multicenter index. Furthermore, in the overwhelming majority of cases these predictions are in full agreement with general expectations.16,33,34 There are only three systems for which significant discrepancies between the expectations and the predictions based on both indices are observed: 23, 32, and 42. In fact, discrepancies between different aromaticity criteria are very common for heteroaromatic species, because in such systems proportionality of the energetic effects of aromatic stabilization and the electron delocalization may depend on the heteroatoms present in the system.12,93–97 Nevertheless, EDDBk index perfectly reproduce the predictions by KMCI, which clearly shows that the bond-orbital projection formalism behind the EDDB approach provides a widely applicable and reliable tool for quantitative evaluation of the multicenter electron sharing and aromatic stabilization effects (at least in the case of organic species). It should be pointed out, however, that the situation might be different in the cases, for which cyclic delocalization of electrons cannot be represented with the resonant covalent forms (e.g. in small charged aromatic rings, metal clusters, etc.) or when the strong cross-ring interactions appear. Then, according to our preliminary (unpublished) studies, EDDBk and KMCI can predict slightly different trends of aromaticity changes; a detailed analysis of this issue will be published elsewhere.
Fig. 8 Summary of the five tests used to compare the performance of EDDBk and KMCI as an aromaticity quantifiers. The complete results with brief comments are collected in Fig. S1–S5 in ESI.† |
The results presented in this work clearly show that the EDDB method safely fulfils the above conditions. For a wide spectrum of organic compounds, the EDDBk index predicts the same trends of aromaticity changes as most of the aromaticity indices from different criteria. What is more, EDDB(r) maps allow for easy identification of resonance-stabilised regions and for tracing of aromaticity changes due to chemical reactions. The comparison of electronic populations from EDDB to other aromaticity measures in the first test proves that the proposed ground-state electron density partitioning in eqn (1) can be considered reliable, and the extracted density layer corresponding to the delocalized bonds (EDDB) in fact seems to be closely related to the properties being characteristic for aromatic rings. In turn, a detailed comparison of EDDB to KMCI shows that for a majority of studied aromatic rings the cyclic delocalization of electrons can be reliably approximated within the framework of the decoupled local resonances at incomparably lower computational cost.
The current implementation of the EDDB method, called RunEDDB, works with both RHF and UHF one-determinant wavefunctions and involves the Hilbert-space partitioning within representation of natural atomic orbitals, which is widely available for most of the popular quantum-chemistry packages through NBO73 and JaNPA101 interfaces. RunEDDB is still under active development and is freely available on the author's website.74
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7cp06114e |
This journal is © the Owner Societies 2017 |