The role of the long-range exchange corrections in the description of electron delocalization in aromatic species

ABSTRACT


Introduction
Quantification of the electron delocalization is crucial for rationalizing electronic structure and molecular properties within the language of such intuitive terms as partial atomic charge, chemical valence, covalency and ionicity, bond multiplicity, etc. [1][2][3][4][5][6] A multitude of different electron delocalization descriptors can be found in the literature, the majority of which are based on the physical-space or the Hilbert-space partitioning of the electron density (ED) from ab initio calculations.In this context, such quantities as the Wiberg's, [7] Mayer's, [8] and the fuzzy-atom bond orders, [9] Bader's diatomic delocalization index, [10] and Fulton's electron-sharing index [11] deserve special attention owing to their reliability and universality. [6,12]Using these descriptors is not limited only to typical two-center bonds.They can easily be extended to cover also the multicenter electron sharing effects (like the three-center two-electron B-H-B bond in diborane) as well as the cyclic delocalization of electrons in aromatic molecules. [13][22] Other popular aromaticity descriptors involving delocalization indices include the para-delocalization index (PDI), [23,24] and the fluctuation index of aromaticity (FLU). [24,25][28][29] In general, the description of electron delocalization depends on the level of the theory used in quantum-chemical calculations. [30]Hartree-Fock (HF) calculations systematically overestimate the charge transfer between bonded atoms due to lack of Coulomb electron correlation in the one-electron density. [10,30]The calculation of delocalization indices at the correlated post-HF levels of theory like the configuration interaction method is rarely performed because it requires correlated density matrices of higher-orders, which are hardly ever available in computational chemistry programs.[32][33][34] However, since the electron-pair density is not strictly defined at the DFT level, the instantaneous Coulomb electron-electron correlation as well as the correlation between electrons of different spin are not taken into account. [35]Consequently, the resulting indices consistently overestimate the electron sharing effects between covalently bonded atoms. [30,36]Furthermore, description of the electronic structure at the DFT level of theory faces yet another obstacle: the incorrect long-range behavior due to the local character of the approximated exchange-correlation (XC) functional. [37]For a long time this issue has been regarded as the principal cause of failure of the KS-DFT approach in terms of the reproduction of such In this paper we address the role of the long-range exchange corrections in description of the cyclic delocalization of electrons in aromatic systems at the density functional theory level.A test set of diversified mono-and polycyclic aromatics is used in benchmark calculations involving various exchange-correlation functionals.A special emphasis is given to the problem of local aromaticity in linear acenes, which has been a subject of long-standing debate in the literature.The presented results indicate that the conventional exchange-correlation functionals significantly overestimate cyclic delocalization of electrons in heteroaromatics and aromatic systems with fused rings, which gives rise to conflicting local aromaticity estimations, especially within the framework of the electronic criteria of aromaticity.In linear acenes, the long-range exchange corrections dramatically improve the performance of almost all aromaticity indices, which after including these corrections consistently predict decreasing local aromaticity going from the inner rings to the outermost ones.molecular properties as barrier heights and reaction enthalpies, [38] excitation energies, [39] van der Waals interactions, [40] nonlinear optical response properties, [41] ionization potentials, [42] etc.45][46][47] Although the local character of the exchange part of the conventional XC functionals has not been found to affect the description of electron delocalization between two covalently bonded atoms, in aromatic rings, in which π-electrons are cyclically delocalized through the net of conjugated bonds over quite large distances, utilizing the LC-XC functionals might be of vital importance. [48]n this paper we address the performance of the LC-XC functionals in accounting for electron delocalization in selected aromatic and heteroaromatic rings as well as polycyclic aromatic hydrocarbons.Since there has been a lot of controversy regarding the problem of local aromaticity in linear acenes ("the anthracene problem"), [49][50][51][52] in this study we take a closer look at the effect of long-range exchange corrections on the description of cyclic delocalization of electrons in these species.

Long-range correction schemes
The long-range correction scheme for XC functionals within the local density approximation (LDA) was originally proposed by Savin. [43]He made use of the standard error function erf to separate the two-electron operator into the short-range (SR) and long-range (LR) parts: where r12 = |r1 -r2|, while µ determines the extent of the short-range region.Applying the standard HF exchange integral to the second term in Eq. ( 1) he got the following definition of the long-range exchange interaction energy: where Ψ iσ represents the ith KS σ-spinorbital.In turn, the LDA exchange functional was combined with the first term in Eq.
(1) giving rise to the following expression for the short-range exchange interaction energy: where In the above equation ρσ(r) is the σ-spin electron density.The corresponding LC scheme for XC functionals from the general gradient approximation (GGA) was developed by Iikura and co-workers. [45]They modified the original Savin's scheme by putting into the averaged relative Fermi momentum included in Eq. ( 4) the gradient terms of the GGA functional.The resulting LC-GGA short-range exchange interaction energy reads Here the GGA exchange-correlation kernel GGA K σ is determined by a particular exchange functional used such that 4/3 GGA 3 ,GGA ( )  .
The ratio between the short-range and long-range terms in Eq. ( 1), µ = 0.33, was determined to reproduce bond distances of homonuclear diatomics up to the third period. [45]ver the last decade increasing attention has been paid to develop the LC schemes for hybrid density functionals, which are the most widely used due to their proven performance for thermochemistry, kinetics, noncovalent interactions, etc. [53,54] One of the most popular LC-XC hybrid is CAM-B3LYP proposed by Yanai and coworkers, who modified the two-electron operator partitioning scheme from Eq. ( 1) using the Coulombattenuating method (CAM). [55]By combining the CAM potential with the B3LYP hybrid [56] they found the optimum portions of the HF exchange to amount to 19% at the short-range and about 65% at long-range (the intermediate region is smoothly described through the standard error function with µ = 0.33).Later Head-Gordon and coworkers came up with their ωB97x hybrid [53] by systematic optimization of Becke's B97 functional. [57]They found that retaining 100% of the Hartree-Fock exchange for long-range electron-electron interactions reduces to a large extent the self-interaction problems associated with global hybrid functionals; at short-range they used the portion of 16% HF exchange with the rangeseparation ratio equals 0.30.In this study the long-range corrected hybrids CAM-B3LYP, ωB97x and ωB97xD (variant with Grimme's D2 dispersion model) [54] are used together with the LC-XC GGA functionals from the LC-scheme by Iikura and coworkers.

Aromaticity indices
The cyclic delocalization of π-electrons gives rise to unusual stability and characteristic structural (reduced bond length alternation) and magnetic (induced ring current in external magnetic field) properties of aromatic species.Thus, apart from the electronic measures of aromaticity based on delocal- (5)  ization indices, in this paper we also investigate the influence of long-range corrections on the aromaticity predictions from structural and magnetic criteria.A brief outline of the descriptors utilized in the present study is given below.Please note that all the provided definitions of electronic measures of aromaticity apply to the closed-shell systems only, although analogous expressions exist for open-shell species.

NICS
The nucleus-independent chemical shift (NICS) [58][59][60] quantifies diatropicity/paratropicity of the induced ring current by means of the effective magnetic shielding measured usually at the centroid of aromatic ring in the external magnetic field, NICS(0), or x angstroms above/below it, NICS(x).The axial component of the nucleus-independent chemical shift 1 angstrom above the ring centroid, NICS(1)zz is one of the most preferable measures of local π-aromaticity since it usually gives aromaticity predictions consistent with results based on other criteria.The more negative value of NICS, the more aromatic is the molecular ring in question.Contrariwise, highly positive values of NICS usually indicate antiaromatic species.

HOMA
The harmonic oscillator model of aromaticity (HOMA) [61][62][63][64] is a normalized measure of the energetic consequences of deviations of bond lengths {Rα} in the n-membered aromatic ring from the corresponding optimum values o { } R α for an idealized aromatic system: where kα stands for the geometry-to-energy conversion factor.The HOMA index varies from close to 1 for highly aromatic rings, through 0 for nonaromatic rings to negative values for antiaromatic rings.

H RCP
The density of total electron energy at the ring critical point (H RCP ) [65,66] is one of the atoms-in-molecule (AIM) parameters that has been proven to serve as a quantitative measure of aromaticity in a large group of molecular rings (e.g. for a set of 33 phenolic rings of varying aromaticity, and a set of 20 quasi-rings formed by intramolecular hydrogen and lithium bonds, the correlation with HOMA and NICS was found to be very tight: R 2 (HOMA)=0.99,R 2 (NICS(1) zz )=0.95).While collecting data for this study, however, we have found that the density of total electron energy at 1 angstrom above the RCP provides even more reliable description of aromaticity in mono-and polycyclic species than the original index.In particular, standard deviations of H RCP (1) are usually smaller by an order of magnitude than H RCP (0); interestingly, similar observations we have made regarding NICS(0) and NICS(1) but in this case differences between the corresponding standard deviations are much less marked -see Tables S2-25 in the electronic supplementary information (ESI).Therefore, in the present study both variants of this descriptor are used.

FLU
The aromatic fluctuation index (FLU) [24,25] is a measure of uniformity of the kekulean mode of electron delocalization and its bonding difference with respect to an idealized aromatic system.It is defined as follows: where the localization and diatomic delocalization indices read occ , ( ) ( , ), ( , ) 4 ( ) ( ), respectively.Here, the summation runs over all (i,j)-pairs of occupied molecular orbitals (MO) and Sij(Aα) represents the corresponding MO-projections integrated within the basin of atom Aα.The κ index in Eq. ( 10) ascertains that the first term is never lower than 1, i.e.
Atomic overlap integrals, Sij(Aα), of Eq. ( 10) were computed within the fuzzy-atom space (FAS) [9] partitioning scheme.By definition, FLU is near zero for strongly aromatic systems and increases as the molecule departs from the reference (idealized) aromatic system.It has been demonstrated elsewhere that the square root of aromaticity fluctuation index usually displays better correlation with classical aromaticity descriptors. [27]However, to make the trends of FLU changes comparable with all the other indices in this study, the FLU -1/2 has been used in most cases of less aromatic rings.

PDI
The para-delocalization index (PDI) [23,24] measures the average effect of electron delocalization between three para-related atoms in a benzenoid unit.For a molecular ring given by the sequence of bonded atoms A 1 , A 2 , …, A 6 PDI reads [ ] The average two-center index (ATI) [16] has the definition very similar to PDI but instead of diatomic delocalization indices it involves the quantum-mechanical bond orders B α,β , In the above equation M k,l stands for the element of the classical Mulliken matrix within the AO-representation, [67] giving rise to the Mayer's bond order, or the element of the charge and bond-order matrix within the basis of natural atomic orbitals (NAO) [68] leading to the Wiberg's bond order.In this study PDI refers by default to the aromaticity index defined within the FAS-partitioning scheme, while PDI(AO) and PDI(NAO) denote the corresponding average two-center indices.

KMCI
[16][17][18] Within the FAS representation it is calculated as follows: thus being a generalization of the diatomic delocalization index from Eq. ( 10).The multicenter delocalization index, congeneric with KMCI but defined within the framework of the Hilbert-space partitioning scheme, reads In the literature the multicenter index obtained from the physical-space ED-partitioning schemes is usually referred to as IR.However, for the sake of clarity, in this paper we use the KMCI label for all three variants of the index, specifying in brackets the particular ED-partitioning formalism used (if not specified, KMCI refers to the FAS-based version by default, just like in the case of PDI).

EDDB
The electron density of delocalized bonds (EDDB) [20][21][22] is a part of the recently proposed method of partitioning the oneelectron density into components representing different levels of electron delocalization [22] .The EDDB formalism makes use of the recently developed bond-orbital projection technique [69][70][71] to approximate the multicenter electron sharing effects by means of the three-atomic local resonances representing conjugations between the adjacent bonds only.Consequently, the EDDB method is computationally far more efficient than the multicenter indices.In the basis of natural atomic orbitals, {χµ(r)}, the electron density of delocalized bonds is defined as where the EDDB metric in the basis of NAOs reads In the above equation Γ is the standard charge and bondorder (CBO) matrix, Cαβ is a matrix of linear-combination coefficients of the appropriately orthogonalized [72] two-center bond orbitals (2cBO), [70] λαβ is a diagonal matrix of the corresponding 2cBO occupation numbers, and εαβ represents a diagonal matrix of the bond-conjugation factors. [21]For ideally localized (two-center) bond A α -A β εαβ is a null matrix, while in the case of delocalized (multicenter) bonds there is at least one diagonal element in εαβ close to 1. Formal definition of the εαβ matrix falls outside the framework of this brief description.It should be noted in passing, however, that the trace of the EDDB metric matrix from Eq. ( 17) determines the population of electrons delocalized through the system of conjugated bonds denoted by Ω.Two variants of the EDDB function are used in this study: EDDB k (r) and EDDB g (r).The first one describes cyclic delocalization of electrons due to the resonance of the kekulean forms; e.g., for a 5-membered ring (5-MR) we have Ω = {(A1-A2),(A2-A3),(A3-A4),(A4-A5), (A5-A1)}.The resulting EDDB k population does not account for the cross-ring electron delocalization and as such can be regarded as a counterpart of KMCI.The EDDB g (r) function, in turn, considers conjugations between all the chemical bonds in a mole-cule, and, as such, it can be used e.g. to evaluate global aromatic stabilization in large molecular systems. [21]

Energetic and structural properties
At first we investigated structural and energetic characteristics of molecules with cyclically delocalized electrons obtained at the KS-DFT theory level using all the benchmarked exchange-correlation functionals.In particular, we calculated bond lengths at equilibrium geometries and compared them to the corresponding experimental values.The complete results with references to experimental data can be found in ESI -see Figure S1 and Table S1; here Figure 1 summarizes the results only for a few selected systems by displaying experimental (black numbers) and the DFT-calculated bond lengths averaged over both groups of functionals, XC (red numbers) and LC-XC (blue numbers).Green numbers in brackets refer to the average bond lengths from calculations involving only the long-range exchange corrected hybrids: CAM-B3LYP, ωB97x and ωB97xD.In turn, numbers inside rings provides the corresponding percentage errors of the calculated bond lengths considering each aromatic ring separately.Next, we compared the KS π-MO energies determined using BLYP and LC-BLYP with the corresponding HF π-MO energies for benzene, borazine, and thiophene.Additionally, differences between the total KS and HF energies, are presented in Table 1.Similar comparison was made for KS π-MOs in pentacene and the resulting orbital energies together with the corresponding isocontours are depicted in Figure 2. Finally, the benchmark calculations of vertical ionization potentials (IPvert) obtained from the DFT-Koopmans approximation [88] was performed for all the aromatic systems, for which the respective experimental data are available. [89]The calculated percentage errors of IPvert for benzene, borazine, thiophene, and pentacene are presented in Figure 3, while the results for other systems are collected in Figure S2 in ESI.

Bond lengths
Although the long-range exchange corrections have been successful in predicting excited state properties and chemical reactivity, so far they have not been reported to be superior to the conventional XC functionals with respect to molecular structure calculations.In fact, such popular GGA functionals as PBE or the hybrid B3LYP one are still regarded as methods of the first choice for geometry optimizations of the most of organic species despite their incorrect asymptotic behavior in the long-range region. [90,91]The results presented in Figure 1 support this fact but only with respect to the monocyclic aromatics: the average errors of the bond lengths calculated using the LC-XC functionals are more than twice as large as the XC-based ones.In practice, the conventional functionals perform better by about 0.01-0.02Å per bond.However, when only the LC-XC GGA hybrids are concerned, the calculated bond lengths are clearly superior to the ones obtained using the pure GGA LC-XC functionals.Here, the differences between bond lengths from both groups of functionals are much less than 0.01 Å per bond for most of the monocyclic aromatic systems.The situation is quite different in the case of polycyclic aromatic systems, as exemplified in Figure 1 using pentacene and picene.In the former the average error of the calculated bond lengths per benzenoid unit strongly depends on the position of the ring: for the terminal rings the LC-XC functionals seem to outdo the conventional ones while in both the inner and the central ring we get the opposite trend.Again, the GGA hybrid LC-XC functionals provide much better predictions of bond lengths especially for the outer rings where the average error is more than two times lower than from the calculations involving standard XC functionals.Furthermore, according to Table S1d in ESI, the total average errors of the calculated bond lengths for the entire molecule are as follow: 0.8% (XC), 1.1% (LC-XC) and 0.7% (hybrid LC-XC).Thus, both groups of exchange-correlation functionals characterizes similar accuracy of the predicted bond lengths in pentacene, with Handy and coworkers' as well as Head-Gordon and coworkers' LC-XC functionals being at a minor advantage over the rest.A thorough analysis of the results for anthracene and tetracene included in Tables S1b and S1c in ESI indicates that the performance of the GGA hybrid LC-XC functionals is only slightly worse than the conventional ones (the differences are at most 0.1-0.3%)while substantially larger errors produce the pure GGA LC-XC functionals.
In contrast to pentacene, in the case of its angular isomer picene average errors of the calculated bond lengths per ring do not depend on the position of the benzenoid unit and to a large extent also on the type of the exchange-correlation functional (the differences are less than 0.3%).Indeed, total  average errors of the calculated bond lengths regarding entire molecule are as follow: 1.4% (XC), 1.4% (LC-XC) and 1.3% (hybrid LC-XC).The corresponding results for phenanthrene and chrysene presented in Tables S1e and S1f in ESI show that the performance of the CAM-B3LYP, ωB97x, and ωB97xD functionals is principally the same as the conventional ones while meaningfully less accurate bond lengths we obtain from the LC scheme for GGA functionals.

Kohn-Sham orbital energies
The first thing that strikes the eye in Table 1 and Figure 2 is that the long-range exchange corrected BLYP functional gives significantly lower energies of the highest-occupied molecular orbitals (HOMO), as well as considerably higher energies of the unoccupied lowest-lying ones, than the original BLYP one.In fact, the LC-BLYP energies are very much like the corresponding HF MO energies, being only slightly lower in most cases.This comes as no surprise, since the long-range corrected functionals reproduce most of the HF exchange and takes into account also the Coulomb electron correlation (and to some extent the molecular non-dynamical correlation introduced by the GGA exchange). [92]Indeed, a comparison of differences between HF and KS total energies with the electron correlation energies from coupled-cluster calculations suggests that it might be the case: utilizing the LC-BLYP functional gives values reasonably close to the reference correlation energies, while the BLYP functional give quite peculiar results.Moreover, in contrast to the LC-BLYP and HF predictions, BLYP gives rise to negative lowest-unoccupied molecular orbital (LUMO) energies, which is somewhat unreasonable.Although the eigenvalues of the Fock operator referring to the LUMO orbitals are commonly regarded as poor estimations of the electron affinities (EA), [93] one should expect at least the same sign of LUMO and the experimental EA; the latter is positive in all three molecules. [89]igure 3 clearly shows that also vertical ionization potentials approximated by the LC-XC HOMO energies (taken with an opposite sign) are closer to the corresponding experimental values, i.e. the relative errors of IP vert are 2-3 times smaller than those predicted by the HOMO energies obtained with the conventional XC functionals.For all the studied molecules, the LC Iikura scheme overestimates IP vert by on average 10-15%; but there are aromatics in which the error is much smaller, like 2% for pyrrole (see Figure S2h in ESI).This slight overshading of the correlation effects by the LC-XC GGA functionals seems to be fixed in the long-range corrected hybrid ωB97xD, which in more than half of the cases underestimates the experimental vertical ionization potentials by only 1-2%.
In fact, all the LC-XC hybrids are characterized by the best performance in terms of reproduction the IPvert for all aromatic systems in this study.The dramatic impact of the longrange exchange corrections on the energetics of the Kohn-Sham orbitals, which now to a certain extent satisfies the DFT-Koopmans theorem ("fictitious constructs but good approximations all the same"), [94] does not automatically entail changes in the patterns of cyclically delocalized electrons.Indeed, a close inspection of the KS-orbital shapes in pentacene depicted in Figure 2 reveals rather subtle changes introduced by the LC scheme (marked by green rectangles).However, it will be shown in the next section that these seemingly insignificant changes are of vital importance for the description of the cyclic delocalization of electrons.

Cyclic delocalization by means of aromaticity indices
As a preliminary, the EDDB g populations were calculated using BLYP and LC-BLYP functionals for all molecules from the test set, and then dissected into different symmetry components to give an account of the total numbers of σ-and π-electrons delocalized through the net of all the chemical bonds.Table 2 collects the absolute and percentage (in brackets) changes of EDDB g and its components, EDDB g (π) and EDDB g (σ), due to the long-range corrections.A full set of the XC/LC-XC benchmark calculation results of all the other aromaticity indices is provided in Tables S2-26 in ESI.Here Table 3 summarizes those results by presenting the percentage changes of mean values (MV) and the corresponding standard deviations (SD, in brackets) of KMCI, PDI, FLU -1/2 , EDDB k , HRCP(1), NICS(1)zz and HOMA, due to the long-range exchange corrections.These quantities are, in our opinion, the most representative indices for different aromaticity criteria and perform generally better than their other variants (significantly smaller values of standard deviations).Additionally, in Figure 4 the results for picene and pentacene were visualized to provide the qualitative framework and illustrate to what extent the LC scheme can determine the picture of local aromaticity in these systems.In the cases where the effect of LC was much below the accuracy of FLU calculation (affected mainly by the numerical integration accuracy and the theory-level inconsistency of parametrization) the average FLU -1/2 values were assumed to be zero and the standard deviations were not calculated.

π-vs σ-delocalization
Even a cursory glance on the results in Table 2 leads to the conclusion that the long-range exchange corrections do not significantly affect cyclic delocalization of electrons in aromatic hydrocarbons.In fact, a small decrease of the EDDB g (σ) components has actually no impact on the relative changes of EDDB g , but it is in line with also very small but systematic reduction of the CC bond lengths, going from benzene (∆RCC= -0.024Å) to cyclooctatetraenyl dication (∆RCC= -0.026Å).In contrast to σ-delocalization, which is responsible for holding the ideally symmetric structure of the ring, [95][96][97][98] no changes of EDDB g (π) are observed.This particularly should not be surprising since the cyclic delocalization of π-electrons is known to be resistant to any change of the bond lengths that preserve the symmetry of the ring. [99]In pyridine and s-triazine the LC scheme slightly adjusts the effects of π-electron delocalization but the magnitude of the resulting EDDB g changes remains small.In turn, more significant reduction of the electron delocalization is observed for borazine, pyrrole, furan, and thiophene as well as for all polycyclic systems.In all of these systems about 80-90% of the total delocalization is attributed to the π-electron effects, which has a direct bearing also on the EDDB g (π) changes due to the long-range correc-   tions, especially in polycyclic aromatic hydrocarbons.For instance, in pentacene the local character of the exchange part of the BLYP functional leads to overestimation of the π-and σdelocalization effects for about 2.8e -and 0.1e -, respectively.Therefore, for the aromatic systems in this study, there is no need to dissect delocalization indices into symmetry components and in practice the total values of the indices reflect the π-electron cyclic delocalization effects with good accuracy.

Monocyclic aromatics
The results of KS-DFT benchmark calculations summarized in Table 3 indicate that for monocyclic species EDDB k , KMCI, PDI, and FLU -1/2 , follow generally the same trends as the total ED-DB g populations from Table 2.In particular, for aromatic hydrocarbons they consistently predict the kekulean mode of electron delocalization to be invariant with respect to any change of the exchange-correlation functional.In heteroaromatic rings the effect of long-range corrections is roughly the same for all four indices except the borazine molecule, for which KMCI and PDI seem to be much more sensitive to the LC scheme than EDDB k and FLU -1/2 .Although it is rather difficult to strictly explain this discrepancy one should bear in mind that in principle EDDB k and FLU are defined by delocalization indices referring to separately-regarded bonds.Consequently, the opposite changes of electron populations on boron and nitrogen atoms caused by LC to a certain extent compensate each other.In contrast, KMCI quantifies the effect of cooperativity of all bonds in a ring and as such is slightly more sensitive to the electron outflow at any atomic center.PDI is even more receptive to the long-range corrections than KMCI despite the fact that both indices are known to be tight correlated for a wide range of 6-MRs. [16]One should realize, however, that PDI is defined by delocalization indices of weak interactions between para-related atoms, which in nature are rather poorly reproduced by the conventional XC functionals. [40]he last three aromaticity indices are not direct measures of cyclic delocalization of π-electrons, which is somewhat reflected by their quite disorderly and inexplicable values.In particular, for all the systems except borazine NICS(1)zz clearly fails to predict the expected reduction of electron delocalization due to corrected long-range behavior of the XC functionals.The same applies to HRCP (1), but in this case, minor extenuation of the index can be observed in pyrrole, furan, and thiophene.The corresponding changes of aromaticity index based on structural criterion are quite difficult to assess.On the one hand, for most of the aromatic rings we get the results that qualitatively resemble predictions based on delocalization indices.It should be mentioned here that since HOMA is determined directly by the bond lengths it is not surprising that the decrease of σ-delocalization effects is better marked in the results (compare with the last column in Table 2).However, there are three cases in which HOMA fails to give expected predictions: cyclopentadienyl anion, cyclooctatetraenyl dication, and furan.This can partially be attributed to the theory-level inconsistencies in the bond length parametrization [63] and the of the choice of appropriate reference aromatic system, which is known to be a major issue in heteroaromatic rings. [64]Moreover, in the case of furan the dramatic decrease of the standard deviation of HOMA (by about 40%) indicates huge discrepancies between bond lengths calculated using various conventional XC functionals, which may significantly affect the impact of the long-range exchange corrections.

Linear and angular acenes
As follows from the analysis of ΔEDDB g collected in Table 2, using the standard XC functionals leads to significant overcounting of the total population of π-electrons delocalized between all the chemical bonds in angularly and linearly fused benzenoid rings.This effect is very clear for anthracene, where the predicted overestimation of the global aromatic stabilization is by about 50% greater than for phenanthrene.
As the number of the fused rings increases, the difference between the LC effects in the linear and angular acenes becomes less pronounced.On the other hand, Table 3 clearly shows that the impact of the long-range exchange corrections on local aromaticity is more difficult to assess as it strongly depends on the position of the particular 6-MR.The only exception is the H RCP (1) index, for which the influence of the LC scheme is too small to be detectable.
In the case of the angular acenes the only discernible regularity in aromaticity reduction due to the LC scheme considers the central ring in phenanthrene and the inner ones in chrysene and picene.Particularly, EDDB k , KMCI, FLU -1/2 and NICS(1) zz consistently predict loss of aromaticity although its magnitude highly depends on the type of descriptor used; e.g. for the central ring in phenanthrene the aromaticity reduction varies from 5% for FLU -1/2 and NICS(1)zz to almost 40% for EDDB k .PDI and HOMA do not follow this trend displaying increase of local aromaticity up to 13% regardless of the ring position (one should realize however that the changes of HOMA might be statistically insignificant due to the standard deviation changes the same magnitude as the index val-It is worth noticing that the EDDB k populations indicate a relatively small reduction of the cyclic delocalization also for benzenoid units in positions identifying localized Clar's sextets, [100][101][102][103] i.e. the terminal rings in all angular acenes plus the central one in picene.This is in contrast to all other aromaticity indices considered in this study. Contrariwise, for the linear acenes all the local aromaticity indices give qualitatively the same estimations of the LC effect: significant reduction of the electron delocalization in terminal rings is accompanied by the enhanced cyclic delocalization in the central ones.Furthermore, most of the indices consistently predict the aromaticity reduction in terminal rings to increase in magnitude with the number of fused benzenoid units.In other words, the larger the linear acene the larger the systematic error caused by the local character of the approximated XC functional.Although the predicted by HOMA increase of cyclic delocalization in inner and central rings seem to be quite excessive (which is probably caused by relatively large average errors of the calculated bond lengths in these rings -see Figure 1), the results for terminal units are nearly quantitatively the same as the corresponding changes of the NICS(1)zz indices.Even the HRCP(1) index gives qualita- tively the same predictions for aromaticity changes in terminal rings in linear acenes but again, the effects of the LCscheme are too small to be regarded as reliable.
Next, let us examine how the changes of cyclic delocalization caused by the long-range exchange corrections translate into the picture of relative aromaticity of benzenoid units in acenes.In Figure 4 the appropriately renormalized mean values (with the corresponding standard deviations) of local aromaticity indices for particular rings in picene and pentacene are displayed.As follows from Figure 4a, all the indices calculated using the standard XC functionals consistently recognize the inner and terminal cycles in picene as the least and the most aromatic ones, respectively, while aromaticity of the central ring ranks somewhere between them.Admittedly, the large deviations of the values of H RCP (1) and HOMA considerably impede the evaluation of relative aromaticity, but, when considering each XC functional separately, they reproduce generally the same cyclic delocalization patterns as other indices.This result can be easily rationalized by Clar's rule, [100- 103] according to which the covalent form with three π-sextets localized in the central and terminal rings contributes most to the overall resonance hybrid.In turn, even a cursory look at Figure 4b reveals that introducing the long-range exchange corrections does not entail any significant changes in the picture of local aromaticity in picene.The only detectable differences consider slightly smaller standard deviations of values (especially in the case of HOMA) and noticeable reduction of aromaticity in the inner rings, which makes the differences between aromaticity of individual cycles more clear cut.A closer study of Tables S15, S16, S19, and S20 in ESI allows one to draw generally the same conclusions with regard to local aromaticity in phenanthrene and chrysene: all the indices predict terminal rings to be the most aromatic ones regardless of the density functional approximation used.
Linear acenes, have no localized π-sextets and according to Clar's rule their benzenoid units are indistinguishable even though they are not symmetry unique.As follows from Figure 4c, in pentacene this qualitative picture of local aromaticity seems to be partially supported by such indices as PDI, FLU -1/2 , H RCP (1), and HOMA, for which the differences in cyclic delocalization of particular rings are much smaller than the corresponding deviations due to approximated character of the exchange-correlation functionals.On the other hand, this result is in contrast not only with the EDDB k and NICS(1) zz predictions that aromatic stabilization systematically increases when going from outermost to the central ring, but also with KMCI advocating for the exactly opposite trend.However, even a first look at the results in Figure 4d reveals that these contradictory predictions vanish for all the aromaticity indices when the LC scheme is involved: the electronic, magnetic and structural criteria of local aromaticity consistently predict decreasing aromatic stabilization effect from central ring to the outer ones.Furthermore, aside from the results by HRCP(1), the resulting order seems to be independent on the type of the long-range exchange corrected density functional used in calculations, since the differences between local aromaticity of particular rings are statistically significant and well detectable.
These findings hold true also for anthracene and tetracene (see Figure 5 as well as Tables S13, S14, S17, and S18 in ESI).Moreover, the effect of the LC scheme on the electronic indices do not depend even on the choice of the electron density partitioning method, what is particularly illustrated in Figure 5 using KMCI and PDI as an example.[106] It should also be noticed that for anthra- cene and tetracene PDI determined using standard XC density functionals gives essentially the same predictions as KMCI, which is somewhat in contrast to pentacene (see Figure 6 in the text, Figure S3 and Tables S21-23 in ESI).
As indicated in previous subsection, experimental bond lengths in linear acenes are slightly better reproduced in the DFT calculations involving the GGA XC density functionals without the LC scheme.Although the absolute average errors are usually worse only by less than 0.01Å when the corresponding LC-XC functionals are used, discrepancies regarding atomic distances in the inner rings of linear acenes are explicitly better marked.To investigate if and to what extent this minor deterioration of the equilibrium geometries affects the performance of KMCI and PDI, the indices were recalculated using density matrices approximated by the GGA XC functionals but at the geometries from the corresponding LC-XC calculations, and vice versa.The complete results considering all the GGA XC functionals separately are collected in Figure S3 in ESI, while the average values are presented here in Figure 6.The results involving the XC and LC-XC density functionals and the corresponding equilibrium geometries (Figures 6a and 6d) assert again, but more clearly, the total effect of the longrange exchange corrections: noticeable reduction of cyclic delocalization in terminal rings, up 33% (KMCI) and 19% (PDI) in pentacene, and simultaneous enhancement of aromaticity in inner rings, up to 76% (KMCI) and 15% (PDI) in anthracene.Figure 6b indicates that, admittedly, the change of geometry itself enhances cyclic delocalization in inner rings, but the reduction of aromaticity in terminal cycles is much smaller.In fact, the effect of geometry is insufficient to reverse the local aromaticity trends predicted by KMCI, but seems to be adequate enough for the corresponding change of cyclic delocalization patterns from PDI predictions (7% reduction for terminal rings in pentacene and only 3% enhancement for the central one in anthracene).As follows from Figure 6c, the effect of the long-range exchange corrections itself (i.e.without reoptimization of molecular geometries) is clearly superior to the effect of geometry as it is strong enough to correct local aromaticity trends predicted by KMCI.In particular, it leads to reduction of cyclic delocalization in terminal rings, up to 20% (KMCI) and 9% (PDI) in pentacene, and simultaneous enhancement of aromaticity in inner rings, up to 64% (KMCI) and 11% (PDI) in anthracene.Thus, although both electronic and geometry effects are responsive to the LC scheme, it is the former that determines the correct order of aromatic rings in linear acenes, and with this regard KMCI seem to be much more sensitive than PDI.
It is worth mentioning that only the local aromaticity desriptors based on the EDDB formalism or on the magnetic criterion of aromaticity (NICS) yield the correct order of the ring aromaticities in linear acenes.These predictions seem to be independent on the choice of the density functional.This observtion is supported by the congeneric maps of the global electron density of delocalized bonds (EDDBg) and the magnetically induced ring-current density (CRD) from the ipso-centric approach, [107] determined using a standard hybrid GGA functional without long-range exchange corrections, as depicted in Figure 7.

Conclusions
A multitude of aromaticity descriptors from the electronic, structural and magnetic criteria have been used to address the role of the long-range exchange corrections in description of the cyclic delocalization of electrons in aromatic systems at the density functional theory level.A test set of diversified mono-and polycyclic aromatics has been used in benchmark DFT calculations involving various exchange-correlation functionals.The presented results clearly show that the wrong long-range behavior of the conventional XC density functionals may lead to significant overestimation of the cyclic delocalization of π-electrons, which in some cases dramatically changes the overall picture of the aromatic stabilization in a molecular system.The essential conclusions of our study can be summarized as follows: • Regarding reproduction of the experimental bond lengths, conventional XC density functionals seem to be superior to their LC counterparts, but only for the monocyclic aromatics; in the angular and linear acenes the LC-XC hybrids like ωB97x give comparable (or even slightly better) results as the non-corrected functionals.• The long-range exchange corrections affect predominantly the cyclic delocalization of π-electrons, which seems to be particularly well accounted for by the LC-XC GGA hybrids, especially ωB97x and ωB97xD, which also quantitatively satisfy the DFT-Koopmans's theorem for all the studied aromatics.][116][117][118][119][120][121] One should bear in mind, however, that most of them were routinely determined at the DFT level of theory using approximated density functionals without the long-range exchange corrections.Thus, in the light of the results from our investigation, these particular local aromaticity predictions can no longer be considered reliable.][124] As far as the structural criterion is concerned, discrepancies between HOMA and delocalization indices were usually attributed to the problems of parametrization, arbitrary choice of the reference system, and even of the definition of HOMA itself. [51]In general, the above-mentioned issues are important and should not be overlooked, as their clarification may lead to deeper understanding of different manifestations of the aromatic stabilization.The results presented in this paper indicate, however, that in the case of linear acenes incompatibilities between local aromaticity descriptors from electronic, magnetic, and structural criteria originate simply from the uncorrected description of electron density delocalization approximated at the DFT theory level.Accordingly, the postulated multidimensional nature of aromaticity [122][123][124][125][126][127][128] seems not to manifest itself at least in this case.
with the reference correlation energies obtained at the CCSD level of theory,

Table 3 .
Summary of Tables S2-26 from ESI: the numbers before brackets stand for the percentage average changes of aromaticity indices due to the long-range exchange corrections, %∆LC = (MVLC-XC − MVXC )/MVXC×100% (MV -mean value), while the numbers in brackets represent the percentage changes of the corresponding standard deviations.Method: KS-DFT/cc-pVTZ, equilibrium geometries.

Figure 4 .
Figure 4.A comparison of mean values (with the corresponding standard deviations) of selected local aromaticity indices calculated using XC (left column) and LC-XC functionals (right column) for picene and pentacene.The values were multiplied by the power of 10 to be in the range of 0 to 1. Method: KS-DFT/cc-pVTZ, equilibrium geometries.

Figure 5 .
Figure 5.A comparison of mean values (with the corresponding standard deviations) of selected electronic indices of aromaticity calculated using XC and LC-XC functionals for anthracene and tetracene.The values of all indices were multiplied by the power of 10 to be in the range of 0 to 1. Method: KS-DFT/cc-pVTZ, equilibrium geometries.

Figure 6 .Figure 7 .
Figure 6.A Separation of structural and electron-density effects of the long-range exchange corrections on local aromaticity in linear acenes by means of the average KMCI and PDI values (multiplied by 10 4 ).Method: KS-DFT/cc-pVTZ.

•
In the monocyclic aromatic systems all the delocalization indices, except for H RCP , consistently predict the LCscheme related decrease of aromaticity up to about 20% for the molecular rings with heteroatoms; the aromaticity indices based on the magnetic and structural criteria do not follow this trend and their response to the long-range corrections is quite disordered and inexplicable.•In the angular acenes the wrong long-range behavior of the conventional XC density functionals lead to noticeable overestimation of the cyclic delocalization mainly in the central or inner rings.The differences between particular cycles, however, are clear cut enough that in total the LC scheme does not introduce any qualitative changes to the picture of local aromaticity obtained with the standard XC functionals.• In the linear acenes, the local character of the approximated exchange-correlation functionals entails significant overestimation of electron delocalization in the terminal rings, and underestimation in the central/inner ones.It leads to conflicting local aromaticity estimations even within the framework of the same aromaticity criterion (Figure4c); the long-range exchange corrections significantly change the performance of almost all the aromaticity indices, which now consistently predict a decrease of the local aromaticity when going from the inner rings to the outermost ones.
The research was supported in part by the Foundation for Polish Science (FNP START 2015, stipend 103.2015,DS), National Science Centre, Poland (NCN SONATA, grant 2015/17/D/ST4/ 00558, DS) as well as the PL-Grid Infrastructure of the Academic Computer Centre CYFRONET, with the calculations performed on cluster platforms Zeus and Prometheus.MS thanks for the support of the Ministerio de Economía y Competitividad of Spain (Project CTQ2014-54306-P), Generalitat de Catalunya (project number 2014SGR931, Xarxa de Referència en Química Teòrica i Computacional, and ICREA Academia prize), and European Fund for Regional Development (FEDER grant UNGI10-4E-801).DS extends his special thanks to MA, MS, TMK and BP for their mental and financial support as well as constant encouragement without which this work would have never been accomplished.
Experimental and calculated average bond lengths with the corresponding percentage errors per ring (numbers inside cycles); results for each XC separately and references to experimental data are included in ESI.Method: KS-DFT/cc-pVTZ, equilibrium geometries.