How do the Hückel and Baird Rules Fade away in Annulenes?

Two of the most popular rules to characterize the aromaticity of molecules are those due to Hückel and Baird, which govern the aromaticity of singlet and triplet states. In this work, we study how these rules fade away as the ring structure increases and an optimal overlap between p orbitals is no longer possible due to geometrical restrictions. To this end, we study the lowest-lying singlet and triplet states of neutral annulenes with an even number of carbon atoms between four and eighteen. First of all, we analyze these rules from the Hückel molecular orbital method and, afterwards, we perform a geometry optimization of the annulenes with several density functional approximations in order to analyze the effect that the distortions from planarity produce on the aromaticity of annulenes. Finally, we analyze the performance of three density functional approximations that employ different percentages of Hartree-Fock exchange (B3LYP, CAM-B3LYP and M06-2X) and Hartree-Fock. Our results reveal that functionals with a low percentage of Hartree-Fock exchange at long ranges suffer from severe delocalization errors that result in wrong geometrical structures and the overestimation of the aromatic character of annulenes.

1. Introduction

Aromaticity is one of the most important concepts in chemistry, and it is associated with cyclic electron delocalization (or conjugation) in closed circuits giving rise to bond-length equalization, energy stabilization, large magnetic anisotropies and abnormal chemical shifts, among other well-known effects [1,2,3]. Despite the multidimensional character of aromaticity and the recent proliferation of aromatic compounds that extend beyond the organic chemistry realm [4,5], several simple but predictive models have been designed to characterize aromaticity [6,7,8,9,10]. From these models, rules of aromaticity [11] have been obtained that are commonly used in assessing the aromatic character of molecules. Among these rules, the Hückel rule was the first and remains the most widely employed [12,13,14,15].

According to the Hückel rule, annulenes with electrons (n being an integer number) are more stable than the corresponding open-chain polyenes and, therefore, considered aromatic. Although this prediction is met for the first elements of the annulenes series, the rule is soon broken due to out-of-plane geometrical distortions that are more energetically favorable but disrupt the optimal overlaps between p orbitals that give rise to conjugated circuits [16]. These geometrical distortions cannot be considered by the Hückel rule which, inevitably, overestimates the aromaticity of annulenes that are described as planar structures having bond equalization. Conversely, the Hückel rule also predicts that electrons annulenes are unstable with respect to their open-chain analogues, exhibit bond length alternation, and are considered antiaromatic. The application of this rule has not been limited to annulenes and it is currently used in all sorts of molecules, including non-organic molecules [4,5].

Another important rule of aromaticity, complementary to Hückel’s, is the Baird rule [17]. While the Hückel rule is limited to singlet ground states, the Baird rule applies to the lowest-lying triplet and excited singlet states and predicts () -electron molecules to be aromatic (antiaromatic). Aihara calculated the resonance energies of singlet and triplet annulenes, concluding that Baird’s rule was satisfied, i.e., that the lowest-lying excited state of annulenes presented the opposite aromatic character to that in the ground state [18]. Baird’s rule has been called the photochemistry analogue of Hückel’s rule [19] because it has been used to tailor molecules with enhanced photochemical activity [20,21,22].

The most (anti)aromatic molecules do not usually have many atoms in the ring. The optimal overlap between orbitals that favors electron conjugation calls for a particular arrangement of the atoms in the ring that cannot be sustained as the ring size increases. Therefore, both aromaticity and antiaromaticity are expected to decrease with the ring size. Although there is evidence of large aromatic nanorings [23], studies on annulenes have led to the conclusion that systems with more than 30 electrons are nonaromatic [24].

In his seminal paper, Baird studied the velocity at which the antiaromaticity (measured through Dewar resonance energy, DRE) diminished with increasing ring size [17]. Based on a few examples, he found that triplets lost their antiaromaticity with the ring size as fast as singlets [17]. To the best of our knowledge, thus far, there have been very few attempts to actually quantify how the aromaticity rules fade away with the ring size and they were limited to annulenes [25,26]. This is the goal of this paper. We will first perform an analysis of how these rules are expected to change with the ring size as predicted by the model from which they originated. Second, we will study how the geometry distorts these molecules from the optimal (planar) geometry assumed in the model and how it affects the aromaticity. Finally, we will analyze the effect of the computational method in the description of these systems by taking into account several density functional approximations (DFAs) with a varying percentage of Hartree Fock (HF) exchange.

2. Methodology

2.1. Aromaticity Indices

In this section, we will review the expressions of several aromaticity indices based on electron delocalization [2]. In the following, we will assume a molecule having at least one ring structure, which consists of N atoms represented by the string = {,,...,}, whose elements are ordered according to the connectivity of the atoms in the ring. For convenience, we adopt and .

2.1.1. The Aromatic Fluctuation Index: FLU

The FLU index measures aromaticity by comparison with the cyclic electron delocalization of some reference aromatic molecules [27]. Its expression depends on the delocalization index (DI) [28,29,30,31], , which measures the electron sharing between atoms A and B, where is the atomic valence for a closed-shell system [32], and is a simple function to ensure that the first term in Equation (1) is always greater or equal to 1,

is the DI corresponding to an aromatic molecule which has the pattern of bonding . For example, the aromatic reference for C-C bonds is benzene. FLU is close to zero for aromatic species, and greater than zero for non-aromatic or antiaromatic species.

Aromaticity indices based on references, such as FLU or HOMA [33,34], do not measure aromaticity but the similarity with respect to some aromatic molecule. Therefore, they are not adequate to describe reactivity [35,36].

2.1.2. The Bond-Length and Bond-Order Alternation Indices

Two popular indicators of aromaticity are the bond-length (BLA) and the bond-order alternation (BOA), which compare the average of bond lengths and bond orders, respectively, of consecutive bonds in a ring: where and , being the floor function of y, that returns the largest integer less than or equal to y, and x is either the bond length or the bond order. Unlike the atomic charges [37], the bond orders are much less dependent on the computational method and the basis set used [38] and, therefore, aromaticity indices based on bond orders [39] are not highly dependent on the level of theory employed. An exception to this rule is the delocalization error that is present in some DFAs [40] with a low percentage of Hartree-Fock exchange, which causes an overestimation of the aromaticity of some compounds (see Section 3.4).

The definition of Equation (3) presents a serious drawback: it is not well defined for a ring of an odd number of members because its value depends on the order of the atoms in the ring. For the sake of generality, in this paper we adopt an alternative definition of BLA/BOA: and we choose the delocalization index [28,29,30,31] as a measure of bond order.

2.1.3. A Many-Center Electron Delocalization Index:

Giambiagi and coworkers suggested to use the multicenter index, which was previously defined by them to account for the simultaneous electron sharing among various centers [41], as a measure of aromaticity [42]. This index was named and its formulation for single-determinant wavefunctions reads as follows: where is the overlap of molecular orbitals i and j in the atom A, i.e., will provide large values for aromatic molecules. Although it will not be considered in this paper, it is worth to mention that Bultinck and co-workers generalized considering also the delocalization of a non-Kekule arrangement of the atoms in the ring; the index is known as MCI [43]. Some of us have shown that both and MCI are ring-size dependent [44] and, therefore, for convenience, we will calculate the multicenter electron delocalization per atom that can be obtained as .

2.1.4. AV1245 and

Both and MCI are among the less fallible aromaticity indices available in the literature [2,36,45] and, therefore, they have been used in a plethora of cases involving a difficult assessment of aromaticity [46,47,48,49,50,51,52,53,54]. However, these indices present some drawbacks that prevent their use in large rings [55] and, therefore, we have recently designed [55] and tested [40,56] a new electronic aromaticity index, AV1245, based on MCI but free of the shortcomings of this index. AV1245 is defined as the average value of the four-atom MCI index between relative positions 1–2 and 4–5 constructed from each five-atom fragment (five consecutive atoms) along the perimeter of the ring [55]. The latter gives an average picture of the electron delocalization among the atoms in the ring. However, for the purpose of measuring aromaticity, it has been found more adequate to use the minimum absolute MCI value evaluated along the perimeter, [56]. identifies the weakest link, i.e., the fragment with the lowest electron delocalization which is usually responsible for the loss of aromaticity in the ring. A detailed study of all the MCI values along the perimeter has been recently shown to be useful in identifying electron delocalization patterns and it will be the topic of discussion in a forthcoming work. Aromatic molecules are thus identified by large values of ( ) and non-aromatic molecules exhibit very low values. Antiaromatic molecules usually exhibit intermediate values and are more difficult to identify. To this end, the analysis is complemented by the examination of either four-atom MCI profiles along the perimeter of the ring or the BOA, which help differentiate between aromatic and antiaromatic compounds.

2.2. Hückel Molecular Orbital Method

Despite the drastic approximations inherent in the Hückel Molecular Orbital (HMO) approach [12,13,14,15], organic aromatic molecules are usually well described within the HMO method. It is thus usual to learn the HMO method at the same time as Hückel’s rule and other aromaticity measures given by the HMO method, such as the resonance energy (RE), the RE per electron (REPE) or the topological REPE (TREPE) [57]. Some recent works have studied Hückel’s rule from the perspective of electron delocalization [58,59]. Since the studies of Hückel on organic molecules, the concept of aromaticity has extended importantly including all sorts of new aromatic molecules such as metalloaromatic molecules [4,51,52,60,61], fullerens [62], nanotubes [63], porphyrins [64,65] with a Möbius-like structure [66,67,68,69,70,71,72,73], and all-metal clusters [4,5], among others.

For a cyclic polyene of n carbon atoms and N electrons, the Hückel molecular orbital (HMO) method provides a general formulation for its orbitals: where is the atomic orbital 2 of the carbon atom and . Excepting for the first () and the last (), these orbitals are complex and degenerate by pairs (). From these equations, one can easily calculate the atomic overlap matrices (AOMs), the bond orders, as well as many aromaticity indices (see Section 2.1), obtaining analytical expressions [74].

The typical way to assess the aromaticity within the framework of HMO theory is through the study of stabilization or resonance energy. Since the computational chemistry study of resonance energies puts severe limitations to analyze more complicated molecules, especially those containing atoms other than carbon, we prefer to study electron delocalization to assess the aromatic character of compounds [2,55]. In particular, we focus here on the study of which, in the case of HMO theory, takes a very simple form (compare to Equation (5)), where is the bond-order between atoms A and B and it is related to the DI: [74].

In order to compare annulenes of different sizes among them, we will study the normalized quantity [44]. The calculation of the bond orders takes a very simple form for singlet annulenes with electrons [74]. The study of other multiplicities or number of electrons is less evident. The HMO theory does not distinguish between alpha and beta electrons (beyond the fact of permitting only one electron of each kind to populate each orbital) and, therefore, the study of triplets is not entirely satisfactory. However, to a reasonable extent, one can analyze the triplets as cases where one electron is promoted from a l orbital to a orbital. If the orbitals obtained from HMO are those of Equation (7), the results are the same regardless the electrons are promoted from l or orbitals (or whether they are promoted to or orbitals). Hence, the study of triplets can be done without ambiguity. However, the study of neutral singlets, which are expected to exhibit localized structures of symmetry , cannot be conducted because both l and orbitals contribute the same amount to all the bond orders in the perimeter of the ring. In order to enforce the appearance of mesomeric structures, we should obtain a different set of orbitals from the HMO theory. Since degenerate orbitals can be freely combined without altering the energy of the system, it is convenient to recombine each pair of degenerate orbitals among themselves to produce this set of orbitals: for . Unlike those in Equation (7), these orbitals are real and produce one of the two mesomeric structures that one can expect in neutral annulenes, depending on whether the last two electrons occupy a l or orbital. Finally, we find that the study of neutral triplets within the HMO theory is far from obvious because both the singlet-open and the triplet states would be actually described by the same orbital occupancies. For the sake of completeness, we have also included this case in our study. To this end, we have employed the original set of orbitals (Equation (7)) and these results are labeled as . (We could have also chosen the second set of orbitals that lead to mesomeric structures. In such case, we would have obtained different values but the same qualitative trend (see Section 3.1)).

3. Results

3.1. Aromaticity from the HMO Method

First of all, we will study Hückel’s and Baird’s rules from the simple HMO theory [12,14]. Although we cannot expect this theory to provide a reliable description of annulenes (especially the large ones), the results we obtain will provide a high bound to the aromaticity/antiaromaticity expected in these species.

All the values of (Equation (8)) are collected in Figure 1. Interestingly, all the curves conform to the same general expression: , where N is the number of ring members. The values of a do not differ significantly among the four cases giving values very close to the theoretical limit for annulenes, , whereas the values of b present large differences among them: −2.1867, 1.0822, −5.1151, and −5.5910 for , singlet, triplet, and singlet structures, respectively. These results are in agreement with the chemical intuition: (i) for very large annulenes both Hückel and Baird rules break and all the species become equally aromatic, () the initial (anti)aromatic character decreases smoothly with the annulene size. We can also see that singlets, which are aromatic, display values above the limit, whereas triplets and singlets, which are expected to be antiaromatic, present values below this limit. On the other hand, annulenes also behave like antiaromatic molecules, which does not conform with the aromaticity expected in triplets and, therefore, these species are not well described from the delocalization measures we can obtain from the HMO theory. Finally, we can compare the velocity with which the (anti)aromaticity decreases according to Hückel’s and Baird’s rule. The values of b show that the velocity at which the aromaticity decreases with the ring size is about five times smaller for singlet annulenes than the corresponding decrease of antiaromaticity for triplets. Interestingly, the latter is very close to the decrease of antiaromaticity found in singlets, in agreement with the study of Baird on cyclic hydrocarbons using DRE [17]. We have also studied the HMO method forcing bond-length alternation by employing two different resonance integrals for each bond type (single and double) [75,76], finding that the antiaromaticity of singlets decreases more rapidly when the annulenes are forced to exhibit bond alternation. Since we cannot compare singlets and triplets, we cannot extract any relevant conclusion about which rule, Hückel’s or Baird’s, is broken more quickly with the ring size. We can, however, conclude from this data that molecules are more resilient to the loss of aromaticity than to the loss of antiaromaticity, as one would expect from purely energetic grounds.

Figure 1

Values of for the annulenes series against the number of C atoms (N) for different singlets and triplets. The species have been divided according to the number of electrons (4n and 4n + 2) and the spin multiplicity (singlet and triplet).

3.2. Geometrical Relaxation

Thus far, we have studied ideally planar annulenes within the HMO method. In this section, we employ quantum chemistry methods to study the geometry of annulenes, which often do not attain a planar conformation in their ground state configuration [77] (see Figure 2). We will employ HF and three DFAs: B3LYP, M06-2X and CAM-B3LYP. These results will be compared against the benchmark data available in the literature. We will split the results into two blocks: and -electron annulenes.

Figure 2

Geometrical structures of the studied annulenes.

3.2.1. Annulenes

First of all, we have studied the lowest-lying structures of singlet and triplet benzene. The ground state of benzene has been largely studied and its characterization does not present a challenge for DFAs. The molecule is identified by all the aromaticity measures as the most aromatic molecule among the studied annulenes. Conversely, the triplet state of benzene has a more elusive character, presenting two conformations which are very close in energy, the quinoidal (Q) and antiquinodial (AQ) one [17]. The former is characterized by two clear double bounds separated by two radical C atoms, whereas the AQ shows two allylic structures merged by two single C-C bonds. The energy difference between the two conformations is below 1 kcal/mol and they are connected by a transition state with an energy barrier of 2.5 kcal/mol [78]. Both conformers are expected to be antiaromatic, displaying an energy destabilization that is responsible for their prominent photochemical reactivity [79]. The lowest-lying triplet state of benzene has been identified by HF, B3LYP, and CAM-B3LYP as AQ-like, whereas the Q structure could not be located with these methods. Conversely, M06-2X provides a Q-like structure. The aromaticity indices display values which are much smaller than those of benzene and the BOA index reveals a bond-order alternation that is characteristic of antiaromatic molecules. values are small, which in conjunction with the oscillating pattern displayed by the dissected index profile (see Supplementary Material), agrees with the antiaromatic character anticipated in this species.

Singlet [10]annulene presents several low-energy isomers, which are difficult to sort energetically from the results of computational calculations. Most ab initio methods find the twist conformation to be the lowest in energy, the exception being B3LYP, which predicts the heart conformation to be the ground state [16]. Our calculations agree with these results. The most delocalized structure, corresponding to a symmetry, lies more than 30kcal/mol above the ground state geometry [16] and, hence, the isomer expected to be the most aromatic it is actually not stable. On the contrary, the twist isomer (which does not exhibit a planar conformation) presents a very modest aromatic/antiaromatic character according to all the aromaticity indices. In fact, the molecule could be easily classified as nonaromatic from the value, which is very close to zero and, as the dissected index profile shows (see Supplementary Material), the twisted bonds are those that prevent an optimal overlap of p orbitals. Conversely, the heart configuration (lying 0.6 and 2.3 kcal/mol above the ground state according to CAM-B3LYP and M06-2X) is an aromatic molecule as stated by all aromaticity indices. To the best of our knowledge, triplet [10]annulene has been much less studied and there are only examples of triplet dicationic [10]annulene structures on the literature [10]annulene Resulting from Baird’s 4n . Angew. Chem. Int. Ed.. 2015 ">20,10]annulenyl dicationic rings. Chem. Eur. J.. 2016 ">21]. The lowest-lying triplet state of [10]annulene displays naphthalene (B3LYP, CAM-B3LYP and M06-2X) and twist (HF) configurations. Following the Baird rule, this molecule should be antiaromatic. The naphthalene conformation is indeed mildly antiaromatic but the twist conformation is rather nonaromatic (see Table 1).

Table 1

Aromaticity indices for annulenes calculated with HF, B3LYP, CAM-B3LYP, and M06-2X and the 6-311G(d,p) basis set.  values were too small for an accurate calculation of  1/N.

Structure	Multiplicity	Functional	FLU	I_ring¯ 1/N	BOA	BLA	I_ring 1/N	|AV_min |
C_6H_6	S	HF	0.000	0.624	0.000	0.000	0.597	10.25
		B3LYP	0.000	0.625	0.000	0.000	0.603	10.72
		CAM-B3LYP	0.000	0.628	0.000	0.000	0.603	10.71
		M06-2X	0.000	0.626	0.000	0.000	0.603	10.73
C_6H_6	T	HF	0.024	0.393	0.246	0.089	0.341	0.39
		B3LYP	0.025	0.399	0.275	0.090	0.353	1.51
		CAM-B3LYP	0.025	0.408	0.276	0.091	0.363	1.20
		M06-2X	0.041	0.467	0.281	0.056	0.380	0.28
C_10H_10 (twist)	S	HF	0.068	0.421	0.728	0.157	0.339	0.00
		B3LYP	0.052	0.480	0.639	0.128	0.322	0.03
		CAM-B3LYP	0.059	0.466	0.677	0.137	0.341	0.03
		M06-2X	0.058	0.463	0.670	0.136	0.325	0.05
C_10H_10 (heart)	S	HF	0.065	0.436	0.712	0.153	0.380	0.01
		B3LYP	0.000	0.610	0.007	0.009	0.579	5.19
		CAM-B3LYP	0.000	0.614	0.009	0.010	0.579	5.13
		M06-2X	0.000	0.611	0.010	0.010	0.579	5.11
C_10H_10 (naphthalene)	T	HF	0.030	0.470	0.364	0.085	0.353	0.14
		B3LYP	0.023	0.531	0.328	0.064	0.460	0.91
		CAM-B3LYP	0.027	0.524	0.364	0.072	0.446	0.68
		M06-2X	0.028	0.522	0.367	0.073	0.444	0.63
C_10H_10 (twist)	T	HF	0.020	0.478	0.266	0.068	0.295	0.05
		B3LYP	0.020	0.517	0.305	0.066	0.352	0.08
		CAM-B3LYP	0.022	0.513	0.324	0.071	0.347	0.07
		M06-2X	0.022	0.511	0.328	0.072	0.337	0.10
C_14H_14	S	HF	0.050	0.497	0.626	0.136	- ^a	0.49
		B3LYP	0.001	0.605	0.010	0.008	- ^a	4.24
		CAM-B3LYP	0.026	0.561	0.449	0.091	- ^a	1.80
		M06-2X	0.025	0.561	0.438	0.088	- ^a	1.89
C_14H_14 (TS)	S	CAM-B3LYP	0.000	0.609	0.007	0.007	-	4.29
		M06-2X	0.001	0.606	0.007	0.007	- ^a	4.27
C_14H_14	T	HF	0.021	0.502	0.282	0.069	- ^a	0.01
		B3LYP	0.017	0.554	0.302	0.060	- ^a	0.13
		CAM-B3LYP	0.023	0.544	0.348	0.071	- ^a	0.08
		M06-2X	0.022	0.544	0.349	0.071	- ^a	0.04
C_18H_18	S	HF	0.049	0.504	0.616	0.133	0.472	0.57
		B3LYP	0.001	0.606	0.026	0.011	0.573	4.27
		CAM-B3LYP	0.026	0.563	0.446	0.090	0.530	1.80
		M06-2X	0.025	0.561	0.444	0.089	0.529	1.81
C_18H_18 (TS)	S	CAM-B3LYP	0.001	0.609	0.022	0.010	0.572	4.29
		M06-2X	0.001	0.607	0.022	0.010	0.572	4.28
C_18H_18	T	HF	0.018	0.514	0.257	0.060	0.411	0.12
		B3LYP	0.013	0.570	0.265	0.053	0.533	0.58
		CAM-B3LYP	0.019	0.559	0.324	0.066	0.513	0.30
		M06-2X	0.019	0.559	0.324	0.065	0.517	0.46

The ground-state [14]annulene geometry corresponds to a distorted pyrene perimeter [16] that can undergo a facile isomerization reaction through a Möbius antiraromatic transition state [14] annulene requires Möbius antiaromatic bond shifting. J. Am. Chem. Soc.. 2007 ">80]. At the B3LYP level of theory, we locate a minimal structure with symmetry that is aromatic according to all criteria of aromaticity. However, the CAM-B3LYP and M06-2X configurations correspond to a transition state (TS) that through a small energy barrier (ca. 2.5 kcal/mol) connects two -symmetry antiaromatic energy minima with complementary bond-length alternation. The triplet conformer leads to an energy minimum that is only mildly antiaromatic at all the levels of theory.

The Hückel rule predicts that [18]annulene is an aromatic molecule but the situation is similar to [14]annulene. B3LYP finds the molecule to be a highly symmetric () structure that is also aromatic from all aromaticity indices. However, both CAM-B3LYP and M06-2X predict that this structure is actually a transition state that connects two identical structures that exhibit bond-order alternation and, therefore, are antiaromatic. This discrepancy has been attributed to the delocalization error of DFAs with a low percentage of Hartree-Fock exchange, a fact that has been confirmed through the comparison of experimental proton chemical shifts [81]. A detailed study of the delocalization error is done in Section 3.4 of this manuscript. The lowest-lying triplet displays a symmetry (HF) or a symmetry (CAM-B3LYP, B3LYP, and M06-2X). DFAs agree on the triplet [18]annulene being weakly antiaromatic, whereas HF predicts it to be rather nonaromatic.

3.2.2. Annulenes

The photochemical formation of annulenes is very important in excited-state aromaticity [19]. Cyclobutadiene is often used as the paradigmatic example of an antiaromatic molecule following the rule. Although (anti)aromaticity can be easily overestimated with a single-determinant wavefunction [82,83,84], cyclobutadiene is the molecule with the largest bond-length and bond-order alternation and all indicators of aromaticity clearly confirm its antiaromatic character. Baird was first to suggest that triplet -electron annulenes should be regarded as aromatic and confirm it through DRE calculations [17]. In particular, the DRE confirmed that triplet cyclobutadiene is an aromatic molecule. A fact that is further substantiated by its symmetric geometry and the values of the aromaticity indices gathered in Table 2.

Table 2

Aromaticity indices for the studied annulenes calculated with HF, B3LYP, CAM-B3LYP, and M06-2X and the 6-311G(d,p) basis set.  values were too small for an accurate calculation of  1/N.  cannot be calculated for rings with less than six members.

Structure	Multiplicity	Method	FLU	I_ring¯ 1/N	BOA	BLA	I_ring 1/N	|AV_min |
C_4H_4	S	HF	0.101	0.391	0.888	0.249	0.262	- ^b
		B3LYP	0.104	0.416	0.900	0.247	0.266	- ^b
		CAM-B3LYP	0.103	0.398	0.898	0.245	0.264	- ^b
		M06-2X	0.103	0.405	0.898	0.242	0.268	- ^b
C_4H_4	T	HF	0.010	0.507	0.000	0.000	0.433	- ^b
		B3LYP	0.012	0.499	0.000	0.000	0.440	- ^b
		CAM-B3LYP	0.011	0.504	0.000	0.000	0.439	- ^b
		M06-2X	0.011	0.505	0.000	0.000	0.438	- ^b
C_8H_8	S	HF	0.067	0.436	0.726	0.156	0.406	0.30
		B3LYP	0.056	0.477	0.664	0.134	0.441	0.72
		CAM-B3LYP	0.061	0.468	0.693	0.140	0.428	0.52
		M06-2X	0.062	0.460	0.694	0.141	0.427	0.51
C_8H_8	T	HF	0.001	0.590	0.000	0.000	0.534	4.07
		B3LYP	0.001	0.589	0.000	0.000	0.540	4.31
		CAM-B3LYP	0.001	0.593	0.000	0.000	0.539	4.29
		M06-2X	0.001	0.591	0.000	0.000	0.539	4.29
C_12H_12	S	HF	0.063	0.445	0.698	0.153	- ^a	0.04
		B3LYP	0.042	0.511	0.565	0.115	- ^a	0.01
		CAM-B3LYP	0.050	0.494	0.624	0.128	- ^a	0.02
		M06-2X	0.050	0.488	0.624	0.128	- ^a	0.06
C_12H_12	T	HF	0.021	0.487	0.280	0.067	- ^a	0.07
		B3LYP	0.002	0.590	0.033	0.012	- ^a	0.07
		CAM-B3LYP	0.015	0.562	0.288	0.056	- ^a	0.21
		M06-2X	0.008	0.577	0.208	0.039	- ^a	0.13
C_16H_16 (S_4)	S	HF	0.054	0.486	0.651	0.139	0.440	0.33
		B3LYP	0.029	0.551	0.476	0.095	0.513	0.96
		CAM-B3LYP	0.041	0.529	0.564	0.113	0.484	0.64
		M06-2X	0.040	0.526	0.562	0.113	0.484	0.63
C_16H_16 (C_1)	S	HF	0.053	0.488	0.643	0.139	0.452	0.25
		B3LYP	0.029	0.548	0.474	0.096	0.512	0.78
		CAM-B3LYP	0.040	0.530	0.555	0.113	0.487	0.55
		M06-2X	0.039	0.526	0.553	0.112	0.487	0.58
C_16H_16 (C_s)	T	HF	0.001	0.598	0.013	0.005	- ^a	1.08
		B3LYP	0.002	0.596	0.057	0.013	- ^a	1.16
		CAM-B3LYP	0.015	0.568	0.294	0.059	- ^a	0.35
		M06-2X	0.011	0.576	0.252	0.050	- ^a	0.74

Unlike cyclobutadiene, cycloocta-1,3,5,7-tetraene (COT) is not a planar molecule in its ground state. COT shows a boat-like geometry that, although presents bond-length and bond-order alternation, is not a very antiaromatic molecule [16,85]. The aromaticity indices reveal that this molecule could be classified as mildly antiaromatic. Conversely, the lowest-lying triplet state of COT presents symmetry and values of the aromaticity indices that confirm the aromatic character anticipated by the Baird rule. Interestingly, the realization of planar triplet COT in some substituted annulenes has been studied as an acceleration path for the photochemical inversion of the ring [22].

In contrast to smaller annulenes, [12]annulene presents several energy minima in the potential energy surface, five of which lie within 5 kcal/mol according to CCSD(T) calculations [86]. The instability and the easy isomerization [87] of this species explain the difficulty in characterizing it experimentally. All DFAs and HF predict a CTCTCT () geometry with a large bond-order alternation, where C and T stand for the arrangement of C-C bonds that can be either cis (C) or trans (T). Nevertheless, the molecule presents very small values that indicate a nonaromatic character. The lowest-lying triplet state of [12]annulene shows a very similar geometry in agreement with a CTCTCT arrangement of the atoms but corresponding to a CTCTCT () geometry, according to all DFAs. The HF lowest-lying triplet also presents a CTCTCT geometry () but it is severely distorted with respect to the latter one. Regardless of the method used, value is small, prompting us to also classify the triplet state as a nonaromatic species.

As in the previous case, a large number of [16]annulene isomers have been found [88] and, therefore, this species also easily undergoes isomerization [87] even through quantum mechanical tunneling [89]. According to all the methods, the ground state structure can be classified as CTCTCTCT (). Interestingly, there is another structure, which presents a distorted CTCTTCTT structure () lying 5-8 kcal/mol above the ground state structure. The lowest-lying triplet presents a conformation very close to the latter one. Although the singlet is less antiaromatic than COT or cyclobutadiene, we find that both and present a similar character and are more antiaromatic than [12]annulene. Therefore, [16]annulene could be classified as weakly antiaromatic. Likewise, the triplet provides values of the indices that indicate a more aromatic character than [12]annulene triplet. Although values oscillate between 0.35 and 1.16, these values indicate a mild aromatic character in between triplet [12]annulene and triplet cyclobutadiene, which agrees with the fact that this molecule displays a more symmetric structure () than its singlet counterpart ().

3.3. Aromaticity from DFAs

The HMO cannot take into account the geometrical relaxation from a planar structure (permitting an optimal overlap of p orbitals) and, therefore, the results obtained in the Section 3.1 are upper bounds to aromaticity and antiaromaticity in annulenes. In this section, we will consider the optimization of several singlet and triplet annulenes and how it affects the aromaticity of these compounds. As an improvement over HMO results, we will calculate an approximate version of that consists in using only the two-center bond orders of bonded atoms to estimate the N-center delocalization. In particular, we will employ the equivalent of Equation (8) in the framework of quantum mechanics. To this end, we will take and subtract 1.00 to (approximately) remove the sigma contribution to the C-C bond, and multiply all resulting DIs of the bonds around the perimeter of the ring, where we have taken the square root to account for the fact that at the HMO level the DI corresponds to the square of the Hückel bond-order [74]. One can consider the latter as a rough approximation that bridges the HMO definition (Equation (8)) and the actual value (Equation (5)). We have collected the results of in Figure 3.

Figure 3

Values of for the annulenes series in terms of the number of C atoms (N) for the lowest-lying singlets and triplets. The species have been divided according to the number of electrons (4n and 4) and the spin multiplicity (singlet and triplet). Calculations were performed with CAM-B3LYP/6-311G(d,p).

The results show some differences with respect to the HMO values we have shown earlier. singlet compounds show a similar trend to the HMO results except for [10]annulene, which exhibits a smaller aromaticity than expected. This molecule undergoes important geometrical distortions that disrupt the overlap of p orbitals and causes the apparent loss of aromaticity. singlets are expected to be antiaromatic, as is confirmed by the large BOA values (see Table 1 and Table 2), and they become less and less antiaromatic as they increase the ring size (see the increasing values of Figure 3). triplet annulenes follow a very similar trend in agreement with Baird’s rule. Finally, we examine the triplet annulenes that are expected to be aromatic and for which there is no clear trend. On one hand, COT, which displays a planar structure, is actually among the most aromatic compounds, even more aromatic than cyclobutadiene in its lowest-lying triplet state, which actually shows the smallest value among the molecules that are expected to be aromatic. On the other hand, neither [12] nor [16]annulene present a planar structure and, therefore, they exhibit larger BOA values, suggesting that these molecules are rather nonaromatic or slightly antiaromatic. Interestingly, all the aromaticity trends are similar to those found for pure HMO calculations for annulenes and singlets with the only mentioned exception of [10]annulene.

Table 1 and Table 2 also include the actual values of but we have not discussed them because they mostly corroborate the trends we have found with and some of the values for large annulenes could not be obtained with enough precision. This is one of the main shortcomings of for large rings [55], the values become so small that they conflict with the numerical precision of the calculation and, consequently, we cannot obtain a reliable normalization ( ) which permits the comparison among rings of different sizes.

In order to remedy this problem, was recently designed [55]. values for the annulenes series are collected in Figure 4. Although the trends are not so clear as in HMO calculations or the approximate account of aromaticity with , the picture we can extract from does not differ entirely from that we obtained from . The four most aromatic molecules are the same according to both indices and each single-triplet pair shows the same order of aromaticity, e.g., singlet benzene is more aromatic than triplet benzene. The only apparent exception to this rule is [10]annulene that has a larger value for the lowest-lying triplet than for the ground state singlet. The latter is due to the fact that the ground state configuration of [10]annulene is found to be nonaromatic (the second lowest-lying structure, the heart configuration, is actually quite aromatic) whereas the triplet configuration is antiaromatic.

Figure 4

Values of for the annulenes series (lowest-lying singlets and triplets) in terms of the number of C atoms (N). The species have been divided according to the number of electrons (4n and 4) and the spin multiplicity (singlet and triplet). Calculations were performed with CAM-B3LYP/6-311G(d,p).

Finally, we find that FLU, BOA, and BLA give a similar qualitative trend to (see Table 1 and Table 2, and also the Supplementary Material) with the exception of large annulenes ([14]annulene and [18]annulene) that are found to be more aromatic in their triplet state than in their ground state. One should bear in mind that these indices are based on pairwise interactions (FLU is also based on fixed ground-state aromatic references) and they are not as reliable as indices based on multicenter calculations [36].

3.4. The Delocalization Error in DFAs

DFAs suffer the so-called delocalization error [90,91,92], which tends to give delocalized electronic structures. This is the exact opposite behavior of HF, which tends to overestimate the localization of the electronic structures. DFAs with a low percentage of HF exchange are more prone to delocalization errors. This effect has been observed in conjugate systems [93] and aromatic systems [40,81,94,95] including singlet annulenes [81]. Recently, Contreras-García, et al. [96] have suggested that the exact result is often bracketed between two extreme situations, the highly localized one (HF) and a highly delocalized one (LDA, which does not include HF exchange). They have used this result to construct error bars for the calculation of some solid-state properties. Likewise, Burke, et al. [97] have defined the sensitivity of DFAs to the density as the energy difference between DFA calculations from LDA and HF densities.

In this section, we investigate how the delocalization error affects the aromaticity measures in the studied annulene series through the analysis of calculated with four methods that use a different amount of HF exchange. In our study, LDA is ruled out because it is not reliable for gas-phase calculations [98] and B3LYP is used as the most delocalizing method (includes 19% of HF exchange). As the most localizing method we employ HF (that obviously uses 100% of HF exchange), and, as DFA with a large percentage of HF exchange, we have selected M06-2X (54%) and a range-separation functional [91], CAM-B3LYP, which presents a variable amount of HF exchange that goes from 19% at short ranges to 65% at large interelectronic separations.

values for different methods are collected in Figure 5. Our results confirm that in the large majority of cases B3LYP and HF are giving the most and the less aromatic species, respectively, among the DFA studied. In some cases (see, singlet benzene, triplet cyclobutadiene or triplet COT) there are no large differences among the DFAs because these species do not suffer from the delocalization error. However, in many other cases (see singlet annulenes larger than benzene) B3LYP clearly overestimates electron delocalization, sometimes even favoring a different geometrical arrangement of the atoms in the ring that permits larger electron delocalization (see Section 3.2). As expected, HF tends to overestimate electron localization whereas CAM-B3LYP and M06-2X provide very similar results lying between HF and B3LYP values. There is only two exceptions to this rule: the triplet states of benzene and [16]annulene. Unlike other DFAs and HF, the M06-2X geometry optimization of triplet benzene leads to a non-planar structure that does not correspond to the AQ structure [17] and, therefore, it is not as antiaromatic as one would expect. A careful examination of the potential energy surface of the triplet state of benzene confirms that M06-2X does not identify the AQ conformation as a minimum of energy at this level of theory. In the case of [16]annulene, the exception is that HF value is actually larger than that of CAM-B3LYP or M06-2X because, according to HF, the structure is more planar than predicted by these DFAs.

Figure 5

Values of for the lowest-lying states of the studied annulenes obtained with methods using a different percentage of Hartree Fock (HF) exchange: HF (100%), M06-2X (54%), CAM-B3LYP (19–65%) and B3LYP (19%).

Interestingly, singlet molecules (with the exception of benzene) and singlets tend to suffer the most from electron delocalization errors, whereas triplet molecules barely exhibit differences between B3LYP and other DFAs with higher percentage of HF exchange. The fact that M06-2X and CAM-B3LYP present very similar values for most molecules corroborates that it is actually the long-range part of Hartree-Fock exchange the one which is relevant to decrease the delocalization error [40]. In this sense, we recommend the use of a range-separation functional, such as CAM-B3LYP, to study aromatic and antiaromatic compounds.

4. Materials and Methods

The singlet ground state and the first triplet excited state of all the studied structures have been fully characterized. The optimizations have been performed with the Gaussian16 software package [99] using B3LYP [100,101], CAM-B3LYP [102], M06-2X [103] DFAs, and HF in combination with the 6-311G(d,p) basis set [104]. The harmonic vibrational frequencies were calculated at the corresponding level of theory in order to verify the nature of the stationary points of their potential energy surface (minima or transition states).

The calculation of the electronic aromaticity indices (AV1245, , BOA and FLU) uses a QTAIM atomic partition [105] performed by the AIMAll software [106]. The AOM resulting from the QTAIM partition and the molecular geometries are the input for the in-house ESI-3D program [27,31,107], which provides AV1245 [55], [56], BLA, BOA, DIs [29,30], FLU [27], HOMA [33] (included in the Supplementary Material), [42], and MCI [43] values. The numerical accuracy of the QTAIM calculations has been assessed using two criteria: (i) The integration of the Laplacian of the electron density () within an atomic basin must be close to zero; () the number of electrons in a molecule must be equal to the sum of all the electron populations of the molecule, and also equal to the sum of all the localization indices and half of the delocalization indices in the molecule. For all atomic calculations, integrated absolute values of were always lower than 0.001 a.u. For all molecules, errors in the calculated number of electrons were always lower than 0.01 a.u. From our experience, these errors provide sufficient accuracy for all the indices here calculated except for the normalized multicenter indices of large rings, which require a numerical precision of the AOM well beyond the accuracy that one can obtain with AIMall or any other similar software available in the literature. For the latter and other reasons commented in Ref. [55], multicenter indices (MCI and ) cannot be used in large rings.

We have employed Mathematica [108] to perform all Hückel calculations, fittings and extrapolations presented in this manuscript.

5. Conclusions

In this paper, we have studied how the Hückel and Baird rules fade away in cyclic polyenes. According to pure HMO calculations and the seminal Baird study [17], antiaromatic annulenes lose their antiaromatic character at the same speed as the ring size increases, regardless of their multiplicity. However, a two-resonance parameter HMO method, permitting bond-order alternation, shows that the antiaromaticy of singlets decreases more rapidly with the ring size. The conclusion is far less clear from calculations that consider the geometry relaxation because the potential energy surface of annulenes with more than eight carbon atoms often shows several configurations close in energy that display disparate aromatic characters. Nevertheless, density functional approximations reveal that the rules fade away much more quickly than it would be expected from the HMO method; they just do not follow a smooth trend.

The study of the level of theory employed in the calculation of annulenes reveals a clear tendency of density functional approximations with a low-percentage of HF exchange at long ranges to exhibit delocalization errors that lead to the overestimation of the aromatic character of the molecule and sometimes even to wrong geometries. Molecules with large ring structures are more prone to this kind of errors. These results are in line with previous findings [40,81,94] and suggest caution in choosing an appropriate density functional approximation to study aromatic and antiaromatic molecules. In particular, we recommend the use of a range-separation density functional approximation such as CAM-B3LYP.