Electronic Aromaticity Index for Large Rings

We introduce a new electronic aromaticity index, AV1245, consisting in the average of the 4-center MCI values along the ring that keep a positional relationship of 1,2,4,5. AV1245 measures the extent of transferability of the delocalized electrons between bonds 1-2 and 4-5, which is expected to be large in conjugated circuits and, therefore, in aromatic molecules. A new algorithm for the calculation of MCI for large rings is also introduced and used to produce the data for the calibration of the new aromaticity index. AV1245 does not rely on reference values, does not suffer from large numerical precision errors, and it does not present any limitation on the nature of atoms, the molecular geometry or the level of calculation. It is a size-extensive measure with a small computational cost that grows linearly with the number of ring members. Therefore, it is specially suitable to study the aromaticity of large molecular rings as those occurring in belt-shaped M\"obius structures or porphyrins.


Introduction
Chemistry is essentially an experimental science that evolved through experimentation and it has been built upon a series of empirically proved laws. On the other side, quantum mechanics relies on postulates from which a solid theory has been constructed. Both focus on the study of matter, however, quantum mechanics can anticipate the electronic structure of matter and it could, in principle, replace the laws and models of chemistry by physically sound theories. Notwithstanding, after many years of the advent of quantum chemistry, several chemical concepts with high predictable power still prevail. Most of these concepts have not found (and most likely cannot find) a solid root in the quantum theory because there is no observable behind them. One finds many such concepts in the literature (e.g. , chemical bonding, bond order, ionicity, electron population, agostic bond, etc.) that are still widely used to predict or explain the electronic structure of molecules or reaction mechanisms. 1 One of the most employed terms in literature -and one of the most controversial 2 onesis aromaticity. [3][4][5][6][7] Aromaticity is associated with cyclic electron delocalization in closed circuits that gives rise energy stabilization, bond length equalization, large magnetic anisotropies and abnormal chemical shifts, among other effects. Various of these aromaticity manifestations can be measured by appropriate quantities, the aromaticity indices, that allow for aromaticity scales. As a result, nowadays there is a number of indices available in the literature, often offering disparate results about the aromaticity of certain chemical species. 8 The discovery of new aromatic species 9-11 that extend well beyond the realm of organic chemistry has challenged our understanding of aromaticity and it has put forward the limitations of the existing aromaticity indices to deal with these new chemical creatures. 12 We have contributed to calibrate and test the limits of applicability of these measures by designing a series of tests that aromaticity indices should pass. 13,14 The indices based on multicenter electron sharing have been the most successful ones. These indices are free of arbitrary references, can be applied to organic and inorganic species alike and, their latest versions are free of the ring-size dependency of the original definitions. 15 Despite these benefits, the multicenter indices suffer from a series of problems that prevents its application in large rings as those showing in new interesting problems such as the aromaticity of (expanded 16 ) porphyrins 17 and its relationship with UV absorption spectra or Hückel-to-Möbius topological switches. 18,19 The goal of this paper is introducing a new electronic aromaticity index that can be applied to rings of arbitrary size. The work is organized as follows: first, we will review the expressions for multicenter indices and reveal their strenghts and weaknesses; second, we will describe an algorithm that permits to calculate MCI on larger rings but it is still limited to rings of medium size; then we will provide the expression of a new multicenter index that will be finally analyzed in a series of compounds to assess its performance.

Multicenter Indices
In the following we will indicate the coordinates of the electron using the short-hand notation 1 ≡ ( r 1 , σ 1 ) and d 1 ≡ d r 1 dσ 1 for the derivatives. Unless otherwise indicated, we will assume a wavefunction constructed from a single-determinant, for a more general approach we suggest Ref. 20. We will measure the aromaticity in a ring structure that consists of n atoms, represented by the string A ={A 1 ,A 2 ,...,A n }, whose elements are ordered according to the connectivity of the atoms in the ring.
To our knowledge, Giambiagi was first to develop an index to measure multicenter bonding. 21 The same expression applied to a molecular ring was named I ring and used to account for aromaticity 22 where S i j (A 1 ) is the atomic overlap matrix (AOM) of atom A 1 , and φ i (1) is a molecular spinorbital. I ring depends on the order of the atoms in the string for n > 3.
Bultinck et al. 23 suggested MCI, which is computed by sum-ming up all the possible I ring contributions obtained from the permutations of the atoms in A , where P is an operator that acting on A produces all the n! permutations of its atoms. The MCI is related to the n-center electron sharing index (nc-ESI), 20,24-26 δ (A ), by withN A being the particle operator applied to atom A and N A the average number of electrons in A, The r.h.s. of Eq. 4 indicates that MCI(A ) is a mesure of how the electron distribution of n particules in n atoms A 1 , . . . , A n is skewed from its mean and it is related to the simultaneous electron fluctuation among these atomic basins. Therefore, I ring measures the delocalization along the ring while MCI also accounts for the delocalization across the ring. In the author's opinion I ring should be connected to cyclic electron delocalization, which is the term used by Schleyer in the IUPAC to define aromaticity. 3 On the other hand, MCI measures the global electron delocalization that is obviously not limited to the cyclic electron delocalization along the ring structure. This distinction becomes important for small rings (typically four-and five-membered rings) where cross-contributions of the electron delocalization might be more important. 27 These expressions, Eqs. 1 and 3, suffer from ring-size dependency. Indeed, the overlaps entering Eq. 2 are, in absolute value, usually smaller than zero and, therefore, the larger the ring the smaller I ring and MCI values. A few years ago, we suggested a normalization that not only avoids this problem but it also closely matches the classical topological resonances per π electron (TREPE) 28 of aromatic annulenes and their ions. 15 In this paper, for the sake of simplicity, we just retain the normalization that produces ring-size-independent indices by taking the 1/n power of Eqs. 1 and 3, i.e. , I 1/n ring and MCI 1/n . If the value is negative, we will take the power of the absolute value and multiply the resulting number by minus one.

The atomic partition
Different atomic partitions can be used to compute the overlaps in Eq. 2, however, some partitions produce completely wrong results, such as indicating benzene as the least aromatic species of a large set of six-membered rings. 29 This result is in contrast with other electronic aromaticity indices such as FLU, 30,31 which is less sensitive to the atomic partition 32 (see also Table 1) because covalent bond orders do not vary much among atomic partitions. 33 In this respect, the most suitable atomic partition 29,34,35 explored thus far is Bader's quantum theory of atoms in molecules (QTAIM) partition. 36 Three-dimensional atomic partitions such as QTAIM's carry some errors due to the fact that numerical integrations of the AOM are performed in non-regular three-dimensional atomic basins. In this regard, one could advocate for the use of Hilbert-space partitions that involve analytical integrations for atom-centered basis sets. However, it is well-known that Hilbert-space-based population analyses suffer from spurious results when diffuse functions, lacking a prominent atomic character, are included in the basis sets. The results collected in Table 1 suggest that multicenter indices suffer to a greater extent from basis-set dependency than its two-center counterparts. Both I ring and MCI usually produce bogus results for basis sets containing diffuse functions if Mulliken or Löwdin partitions are used to define the atoms in the molecule. Mulliken results are unusually large while Löwdin partition values are too small. The PDI, 37 which measures the para-related 2c-ESI in six member rings, is also largely affected by the inclusion of diffuse functions, whereas FLU, 30 which is constructed from the 2c-ESIs between bonded atoms, only shows important changes for pyridine calculated with Mulliken's partition. Interestingly, for cc-pVDZ basis set, the agreement between Löwdin and QTAIM partition is quite good for all the indices. It is thus apparent that Hilbert space based partitions can only be used in the absence of diffuse functions. Since many chemical systems exhibit highly delocalized electronic structures that call for the use of diffuse functions, we recommend the use of three-dimensional space partitions for the calculation of multicenter indices.

Computational Cost and Numerical Precision
The operation P(A ) in Eq. 3 produces n! summations and thus constitutes the bottleneck of the MCI calculation. Overall, if carried out naively, the computational cost of MCI goes like m 3 n!/2, where m is the number of occupied orbitals and we have taken into account that there is only (n − 1)!/2 distinct permutations for the elements in A .
The computation of three-dimensional atomic partitions carries numerical precision errors in the calculation of the AOMs and the total MCI error grows with the ring size. There is three effects as ∆n members are included in A : (i) on one hand, the number of summations in Eq. 3 increases by roughly a factor of n ∆n , (ii) each term in the summation of Eq. 1 involves ∆n additional multiplications and (iii) the enlargement of the ring size probably carries the increment of the dimension of the AOM, as larger molecules have more electrons and, therefore, more orbitals will be involved in the calculation. AOM elements are much smaller than zero for the delocalized orbitals that contribute to the MCI and hence, as the ring size increases the error in each summation can easily grow beyond the machine precision. In addition, for calculations using correlated wavefunctions (see next section) this situation worsens considerably. For single-determinant methods (or Kohn-Sham DFT) the dimension of AOM is proportional to N 2 , where N is the number of electrons, but for correlated calculations the dimension goes like m 2 , where m is the number of occupied orbitals.
There is a number of tips we can follow to improve the computational performance and also minimize the error in the MCI. First, we should use very accurate numerical integrations schemes to avoid large error propagations. Second, we should make a careful selection of the orbitals involved in the MCI calculation. Bultinck et al. introduced the pseudo-π method, 38,39 which reduces the size of m by considering only the π orbitals in the system. The approximation works very efficiently for polycyclic aromatic hydrocarbons but many large rings, such as belt-shaped Möbius molecules, are not planar and do not render themselves to a classification of its orbitals into σ and π types. One can also manually select the set of orbitals that can contribute to the aromaticity of the system. Usually only valence orbitals are delocalized enough to contribute to the MCI and, therefore, the cost can be safely reduced by only including these orbitals. For highly symmetric molecules, the use of symmetry constraints can also reduce the cost by considering only symmetrically distinct permutations. Besides, the implementation of Strassen algorithm for matrix multiplication can also contribute to lower the cost to m 2.8 (n − 1)!/2.
Overall, the large error propagation cannot be avoided and it is an inherent problem for the calculation of MCI and I ring ; fortunately, these problems only show for large rings (n ≥ 10).

Correlated wavefunctions
The calculation of the MCI from a correlated wavefunction through Eq. 4 bears a large computational cost, associated to the computation of the n-density. The calculation of n-order densities using correlated wavefunctions for n > 3 is typically beyond the capabilities of current computational resources, except for very small molecules. For this reason, many authors, 40 including ourselves, 15 have decided for approximations to the n-order densities that provide cost-efficient MCI calculations using correlated wavefunctions. Lately, we have introduced two approximations of nc-ESI that can be casted in the following expression for the MCI: where a is a constant and {n i } are natural orbital occupancies. a = 1 corresponds to the single-determinant expression of the ndensities, whereas a = 1/2 and a = 1/3 are our proposals. The latter has been shown to provide the most accurate results for 3c-ESIs, even better than more sophisticated 3-density approximations. 20 However, the exact calculations of 3c-ESIs tend to provide very small values and certain capabilities of the 3c-ESIs are lost upon inclusion of electron correlation. 20 For instance, Eq. 6 with a = 1 is the only variant of MCI calculation that permits to distinguish between 3c-2e and 3c-4e interactions. 24,41 For this reason, in recent aromaticity studies using correlated wavefunctions we have preferred the a = 1 approximation, which provides more sensible numbers. 15,42

A new algorithm for MCI
In order to bring down the computational cost it is much more profitable to reduce the number of distinct permutations in A than reducing the scaling of the number of basis functions. To this aim, we devise an algorithm that screens the superatomic overlap matrices (SAOMs hereafter) that result from the partial summation of some indices, e.g. , These matrices can be precomputed only for the SAOMs that produce some significant (above some threshold) value. The reduction of the number of SAOMs decreases the cost of the MCI calculation by decreasing the number of possible permutations. The most efficient way to construct these matrices is by successive duplication of the SAOM order following the binary representation of the number of atoms in the ring (n). First of all, we construct a supermatrix of order p, where p is the largest number 2 p that is smaller or equal to n, by consecutive self-combination of lower order matrices. Namely, we first construct the SAOMs of order 2 (2-SAOMs) from the AOMs of the atoms in our ring and the 2-SAOMs are combined among themselves to build the 4-SAOMs and so on until 2 p -SAOMs are formed. Finally, the pertinent low-order matrices (generated through this process) are added to 2 p -SAOMs until the final n-SAOMs are constructed. The MCI is calculated by summing up the traces of all the n-SAOMs. Let n=∑ p i=0 a i 2 i be the binomial representation of n. The number of steps in the algorithm equals p + ∑ p−1 i=0 a i . The strengh of the algorithm is that a screening is applied to every step of the procedure so that the number of SAOMs at each step is reduced. The algorithm could be enhanced by the use of an orbital localization scheme to increase the sparsity of the SAOMs in Eq. 6.

A new electronic aromaticity index
Despite the speed up gained by the algorithm described in the last section, we anticipate that the calculation of MCI on non-planar rings of more than 16 members still poses a serious challenge. Moreover, for such large rings the numerical error in MCI and I ring is quite big. For these reasons, we prefer to explore the possibility to use other electronic mesures to account for the aromaticity of large rings. The FLU has been successfully used in the past to study the ground-state structure of several porphyrins. 17 However, we have discussed the limitations of reference-based indices such as FLU to account for the aromaticity in molecules with stretched bonds or to recognize the aromaticity of transition state structures. 13,43 In a previous work 44 we studied the long-range delocalization patterns in annulenes and found that the interactions separated by an even number of atoms (e.g. , 1,4 and 1,6) are larger than those separated by an odd number of atoms (1,3 and 1,5, respectively) in aromatic molecules. This finding is in line with the well-known fact that the meta-interaction is smaller than the para-interaction in benzene 45 (and most aromatic six-membered rings 37 ) despite the larger separation of the atoms in para position. Since we want to capture multicenter delocalization effects, we suggest to average all the 4c-ESI values along the ring that keep a positional relationship of 1,2,4,5, as indicated in Figure 1. This quantity (hereafter named AV1245), will be large if there is an important delocalization between bonds 1-2 and 4-5, a condition that only occurs in conjugated circuits. Notice that the atom in position 3 is skipped intentionally because we want to include two 1-to-4 interactions (1,4 and 2,5) and measure the extent of transferability of the delocalized electrons in the 1-2 bond to bond 4-5 and viceversa. We have explored other possibilities including these interactions (e.g. AV124) but none have given as good results as this one. AV1245 does not rely on reference values, and it does not present any limitation on the nature of atoms, the molecular geometry or the level of calculation (we suggest to use Eq. 4 with a = 1). Besides, it does not suffer from large numerical precision errors, size-extensivity problems or unfavorable computational scaling. AV1245 bears a very modest cost of 12m 3 n that thus only grows linearly with the ring size. The only limitation of AV1245 is that it cannot be applied to rings of less than six members. However, for such small rings we can safely use MCI.

Computational Details
We have studied several molecular systems that can be classified into five sets. The first set contains a series of neutral annulenes: benzene, C 8 H 8 (planar), annulene [10], annulene [12], annulene [14], annulene [16] and annulene [18], which will be used to assess the capabilities of the indices to recognize aromaticity as the number of pi electrons increases. The second set contains a series of molecules that are used to validate the performance of the new index in very different carbon-carbon bonds: cyclohexane, cyclohexene, benzene and the transition state of the Diels-Alder reaction occurring between butadiene and acetylene. A series of heteroaromatic compounds with six-membered rings (6-MR) will be used to guarantee a correct description of molecules containing heteroatoms: pyridine, pyridazine, pyrimidine, pyrazine, triazine, borazine and purine. The fourth set comprises some polycyclic aromatic hydrocarbons: naphthalene, anthracene, phenanthrene, cyclobutadiene, azulene and pentalene for which we analyze the 6-MRs and the peripheral rings. For the sake of completeness we have included two non-aromatic seven and eight-membered rings molecules (C 7 H 14 and C 8 H 16 ), an aromatic seven-membered ring (C 7 H + 7 ) and hexaethynyl-carbo-benzene (C 18 H 6 ) as the molecule with a very large ring. 46 This gives a total of seventeen 6-MRs, three 7-MRs, four 8-MRs, one 9-MR, four 10-MRs and one 12-MR (total 30 points) that will be used to compare the performance of MCI and AV1245. In addition, there is four 14-MRs, one 16-MR and two 18-MRs data for which only AV1245 values will be provided.
The calculations have been performed at the B3LYP/aug-cc-pVDZ level using the g09 package. 47 The AOM calculations use the QTAIM partition and have been computed with AIMall package 48 and inputed into ESI-3D code. 49

Results
In Figure 2 we plot AV1245 against MCI for the 6-MR studied in this work. The data comprises aromatic and non-aromatic rings, the transition state of the Diels-Alder reaction (which is a difficult case for several aromaticity indices) 43 and the set of heteroaromatic molecules. AV1245 gives an excellent linear correlation with zero intercept and a Pearson coefficient of r 2 = 0.98. The only small discrepancy between both indices is pyrazine which according to AV1245 is slightly more aromatic than benzene.
In Figure 3 we plot AV1245 and MCI 1/n against the number of carbon atoms for the series of neutral annulenes (set 1). The number of π electrons equals the number of carbon atoms and, therefore, according to Hückel's rule, one expects a zig-zag plot where aromatic and antiaromatic compounds alternate. Although some antiaromatic compounds present negative MCI values, the experience indicates that many antiaromatic compounds simply give very small number and are barely distinguishable from nonaromatic compounds. 15 Let us now analyze the ring-size dependency of both indices by including all rings from 6-MR to 12-MR in the plot of Figure 4. The indices do not show a linear correlation trend but, instead, a power-law dependency that fits most of the data points studied.  A few exceptions have been marked in blue and red. The blue bullets indicate antiaromatic molecules that show negative values for both MCI and AV1245. One is thus deemed to conclude that, as it happens in MCI, antiaromatic compounds exhibit either negative or very small AV1245 values. The red bullets correspond to molecules with an odd number of ring members. These molecules clearly deviate from the general trend and put forward that the conjugation mechanisms in such rings is less obvious. The aromaticity in these compounds is severly underestimated by AV1245 and until further explorations are performed AV1245 should be used in rings of even number of atoms.
Finally, we collect some AV1245 values for large rings in Table 2. Among the polycyclic aromatic hydrocarbons (PAHs), the most aromatic 6-MR is the external ring of phenanthrene, whose two-ring naphathalenic structure is the least aromatic 10-MR among PAHs. The peripheral ring of naphthalene shows similar AV1245 value to the annulene [10] structure, whereas neither the peripheral ring of anthracene or phenanthrene are as aromatic as annulene [14]. Conversely, hexaethynyl-carbo-benzene also exhibits similar aromaticity to annulene [18], which confirms the prominent aromatic character of this molecule. 46

Conclusions
We have introduced a new electronic aromaticity index, AV1245, consisting in the average of the 4-center MCI values along the ring that keep a positional relationship of 1,2,4,5. AV1245 measures the extent of transferability of the delocalized electrons between bonds 1-2 and 4-5, which is expected to be large in conjugated circuits and, therefore, in aromatic molecules. AV1245 analyzes the aromaticity of individual rings (local aromaticity) but it could also be used to measure the extent of all 1,2,4,5 conjugation interactions in a molecule and, therefore, account for the global aromaticity of annulated species such as benzenoid macrocycles. 50 A new algorithm for the calculation of MCI for medium-sized rings has been introduced and used to produce the data for the calibration of the new aromaticity index.
Our results indicate that AV1245 correlates very well with MCI excepting for rings with an odd number of members. AV1245 does not rely on reference values, does not suffer from large numerical precision errors, and it does not present any limitation on the nature of atoms, the molecular geometry or the level of calculation. It is a size-extensive measure with a small computational cost that grows linearly with the number of ring members. Therefore, it is specially suitable for studying the aromaticity of large molecular rings as those occurring in belt-shaped Möbius structures or porphyrins.