Three Queries about the HOMA Index.

HOMA (Harmonic Oscillator Model of Aromaticity) is a simple, successful, and widely used geometrical aromaticity index. However, HOMA can also be used as a general molecular descriptor appropriate for any type of molecule. It reaches the global maximum for benzene, whereas the potent magnetic aromaticity NICS index has no lower or upper limits. Hence, questions arise and go beyond mere differences between the geometric and magnetic aspects of aromaticity: (1) Does a molecule of aromaticity greater than that of benzene, but undisclosed by the HOMA definition, exist? (2) Can the Kekuléne cyclohexatriene moiety with HOMA = 0 exist as a part of a larger system? (3) Can the geometrical aromaticity index be defined better? Our answer to the first query is "It is not likely enough", to the second, "Why not define HOMA using a less mysterious molecule than cyclohexatriene?", and to the third, "It is possible to construct another fair geometrical index, but is it better for evaluating aromaticity?" To find these answers, we have studied: (1) the HOMA and NICS indices of over 50 hexahomosubstituted benzenes, (2) the HOMA, as well as EN and GEO, indices of over 100 triply fused hexasubstituted benzenes, and (3) the HOMA and new Geometrical Auxiliary Index (GAI) , of different unsaturated and saturated, aromatic and aliphatic hydrocarbons including all alkane constitutional isomers composed of up to nine carbon atoms.

Introduction

Modern quantum-mechanics-based computational chemistry is focused on obtaining exact data. However, for the vast majority of molecules, obtaining exact data is still unachievable. Moreover, for most computational chemists, only approximate theoretical levels are accessible. On the other hand, there are plenty of organic physical chemistry (semi)empirical descriptors, such as substituent-, solvent-, reactivity-, aromaticity-, topological indices, etc., that may be helpful for directing the research and modeling of the properties of interest. Even though it seems that interest in those descriptors already reached its highest intensity in the 1950s–1980s, such indices are helpful in categorizing and preselecting species for further scrupulous quantum-chemical analyses. Despite the fact that they do not always have a clear physical meaning, their use is easier, quicker, cost-effective, and more available to a broader chemical community than the robust computational methods. Such simple parameters can be auxiliary in materials design, and thus they deserve to be discussed, understood, and developed.

HOMA (Harmonic Oscillator Model of Aromaticity) is one of the simplest, most successful and most widely used aromaticity indices.[1,2] It embraces the geometrical aspect of aromaticity in a formula in which benzene bond length is the internal standard of perfect aromaticity. The more different from benzene and the more unequal and alternating the bonds are, the lower the HOMA index of the analyzed ring is, denoting how less geometrically aromatic it is. The HOMA index can be expressed as followswhere R and Ropt stand for the ith bond length in the analyzed ring and the reference benzene ring (1.388 Å), respectively, n is the number of the CC bonds in the ring, and α = 257.7 Å–2 is a normalization factor making the unitless HOMA index equal to 1 for perfectly aromatic benzene and 0 for a perfectly alternating hypothetical Kekulé cyclohexatriene ring.

The HOMA index is parameterized for the most important heteroatoms and can be reparameterized further.[3,4] It can be decomposed into purely geometrical and purely energetical components.[5,6] Moreover, it has recently been shown that HOMA can be used also as a general molecular index adequate for unsaturated or saturated, cyclic or acyclic, linear or branched (hydrocarbon) molecules.[7,8] We proposed calling this generalized property the savoricity. Quite recently we have additionally shown that it can be a topological index in the structural formula version of the graph theory.[9] The HOMA index has also been a source of several interesting modifications.[10,11] The multidimensionality of aromaticity and the incompatibility of different aromaticity measures have been deliberated for a long time.[12−16] However, recently, it seems that direct assessment of aromaticity is approached. For example, the impact of the π-electron cyclic delocalization on aromatics stabilization energy was quantitatively accounted for by determining the atom–atom charge transfer, repulsion energy, and spin orbitals population at π-electron sites.[15]

The unconformity between the indices can be to some extent associated with their different construction. Indeed, HOMA is a relative index reaching its maximum for benzene, whereas the definition of the other most extensively used NICS (Nucleus-Independent Chemical Shift) magnetic aromaticity index[17−21] sets out no lower or upper limits of this index. As a consequence, juxtaposition of the HOMA and NICS values raises questions not only about the differences between the geometric and magnetic aspects of aromaticity. There are numerous molecules with NICS(0) or NICS(1) indicating aromaticity larger than that of benzene, so: (1) Does a molecule of aromaticity greater than that of benzene exist, but is undisclosed by the HOMA definition? A hypothetical cyclohexatriene molecule has by definition HOMA = 0. However, (2) Can the Kekuléne cyclohexatriene moiety with HOMA = 0 exist as a part of a larger system? The third question that appears is about the construction of the geometrical aromaticity index: (3) Can the geometrical aromaticity index be defined better? This paper attempts to discuss these three questions. These questions do not exhaust the issues that might be raised. One of them, raised by one of the reviewers of this paper, is Does the HOMA index can be reformulated to well distinguish nonaromatic from antiaromatic structures? However, the answer to this problem is difficult and it is worth to be a central issue of the future study.

Calculations

All calculations were performed using the B3LYP functional,[22,23] the 6-31G** Pople type basis set,[24] and Gaussian09 software.[25] The 6-31G** basis set was shown to perform fairly well in geometry, frequency, energy, and electron density calculations.[26,27] Each minimum was confirmed by checking that all harmonic frequencies were positive and that the optimized molecules exhibited proper symmetry. The alkane conformers were found using the molecular mechanics search of conformer distribution as implemented in Spartan’14 suite of programs.[28] The conformers were reoptimized at the B3LYP/6-31G** level, and all were local minima on PES.

Results and Discussion

Does an “Overaromatic” Carbon Ring Exist?

It is clear that search for a molecule with the aromaticity exceeding that of benzene cannot be based on the HOMA index alone: it reaches the global maximum for benzene by definition (eq ). Hence here, the NICS(1) or NICS(0) indices serve us as preselection parameters: if they indicate that the magnetic aromaticity is greater than for benzene, then the system is analyzed further. The NICS(1) index is a reverse chemical shift value calculated 1 Å above the center of the ring, while NICS(0) is analogous value taken in the ring center.[17−19] We assume that “overaromatic” organic molecules can be found primarily among the planar six-membered rings with equal C–C bonds. Although this condition seems to be necessary, it is not sufficient. In fact, all C–C bonds in the central coronene ring are equal, but this ring exhibits a decreased aromaticity (HOMA ≈ 0.75) and the rings in (nonplanar) [6]radialene and (planar) cyclohexanehexone (CHONE) are nonaromatic and strongly nonaromatic (HOMA ≈ −1.3 and −5.4, respectively).[29]

Collation of the NICS(1), NICS(0), and HOMA aromaticity indices and the CC bond distances for about 50 hexahomosubstituted benzenes (Table ) contains both existing and hypothetical molecules. A large set of such compounds was not analyzed earlier in terms of their aromaticity. Moreover, although we calculated several other hexahomosubstituted benzenes, we skipped them because we received no significant new information. The hypothetical structures, like the hexalithium-, hexaphosphaethyno-, hexaisocyano-, hexasiloxane-, or hexachlorosylbenzene, are used here to provide smooth changes of the substituent properties. Nevertheless, some of them are probably now synthesizable. Several molecules, e.g., F, O–, CN, or C≡CH hexasubstituted benzenes, exhibit perfect D6 symmetry. But, most of them have a reduced symmetry because of substituent bulkiness, steric repulsion, and interactions which can even lead to losing ring planarity. Nevertheless, all studied compounds exhibit only one ring CC bond length, and even if they have decreased symmetry, they experience no NICS(1) index splitting.[21]

Table 1

NICS(1)ZZ, NICS(1), and NICS(0) (ppm) and HOMA Aromaticity Indices, and d(CC) Bond Distances (Å) in Hexahomosubstituted Benzene Rings Estimated at the B3LYP/6-31G** Levela

subst/mol	NICS(1)zz	NICS(1)	NICS(0)	d(CC)	HOMA	subst/mol	NICS(1)zz	NICS(1)	NICS(0)	d(CC)	HOMA
OK	94.373	27.606	35.049	1.542	–4.457	SMe	–21.390	–10.798	–8.705	1.414	0.918
ONa	48.891	10.847	14.473	1.477	–0.699	OCOMe	–21.541	–10.191	–13.176	1.401	0.994
CHONE	2.699	4.169	16.138	1.554	–5.416	OMe	–21.988	–11.513	–14.818	1.403	0.988
coronene	1.737	–5.362	–1.233	1.428	0.740	OCHO	–21.989	–10.360	–13.491	1.401	0.994
Li(*)	–2.641	–0.258	–31.400	1.413	0.931	OC≡CH	–21.995	–11.700	–15.893	1.397	1.000
OBr	–13.252	–9.800	–10.999	1.414	0.916	BeH(*)	–22.183	–10.208	–3.373	1.434	0.633
CP	–13.288	–8.943	–9.531	1.427	0.756	OCN	–22.215	–11.882	–16.357	1.396	1.000
SH	–16.060	–8.638	–7.735	1.411	0.944	CClO	–22.576	–10.976	–13.070	1.405	0.981
C≡CF	–17.236	–8.968	–8.260	1.416	0.899	CF3	–22.768	–11.206	–11.598	1.410	0.951
OCl	–17.323	–10.527	–12.168	1.409	0.960	OPh	–22.823	–8.727	–11.442	1.402	0.991
O–	–17.327	–9.792	–11.823	1.506	–2.106	F	–23.214	–11.870	–17.214	1.393	0.997
OF	–17.328	–9.127	–11.309	1.420	0.854	CFO	–23.233	–10.668	–11.408	1.404	0.985
C≡CH	–17.428	–9.196	–8.817	1.417	0.889	Ph	–23.304	–8.261	–7.113	1.413	0.931
OSiH3	–18.256	–10.755	–14.933	1.400	0.996	SiH3	–23.510	–10.065	–5.392	1.419	0.866
Br	–18.321	–9.478	–10.000	1.403	0.988	Me	–23.745	–10.380	–8.750	1.409	0.958
OCaF	–18.603	–11.677	–15.531	1.415	0.909	COONa	–23.882	–10.393	–11.947	1.409	0.956
NC	–18.728	–10.350	–13.267	1.410	0.951	COOH	–24.169	–10.835	–11.098	1.412	0.936
OCaH	–18.837	–11.355	–14.849	1.416	0.901	Li(☆)	–24.204	–7.692	–5.128	1.422	0.829
BH2	–18.887	–7.712	–3.782	1.421	0.842	COO-	–24.236	–11.140	–12.182	1.434	0.624
OMgH	–19.662	–11.652	–15.207	1.397	1.000	NH3+	–24.370	–11.429	–14.043	1.431	0.688
Cl	–19.857	–10.127	–11.216	1.404	0.984	COOLi	–24.799	–10.705	–12.523	1.410	0.953
NH2	–20.125	10.279	12.645	1.401	0.993	MgH	–25.034	–9.402	–0.902	1.425	0.787
OH	–20.192	–10.835	–14.983	1.395	1.000	CHO	–26.496	–12.240	–11.804	1.410	0.951
ONO2	–20.244	–8.829	–9.951	1.404	0.984	C6(1)	–27.120	–12.890	–25.968	1.328	–0.199
CN	–20.472	–10.418	–10.857	1.411	0.944	H	–29.191	–11.305	–9.846	1.396	1.000
OLi	–20.613	–11.672	–15.749	1.410	0.951	BeH(☆)	–36.829	–13.614	–4.422	1.410	0.951
SLi	–21.205	–10.194	–10.468	1.421	0.842

Systems are ordered according to the decreasing NICS(1)ZZ value. For some unusual substituents, see Scheme . (*) denots bonding with the C-atoms. (☆) denotes bonding to the C=C bonds, (1) a molecule.

Scheme 1

General Structural Formula of the Hexasubstituted Benzenes and Some Structures Which May Be Incorrectly Recognized Based Only on Table Notation of Substituents

Virtual BeH and Li hexasubstituted D6 benzenes. At the B3LYP/6-31G** level, the (*) types exhibit a few imaginary modes while the (☆) types show all real frequencies.

Among the molecules collected in Table , there are some exceptional cases (Scheme ). The C6 molecule converges to the D3h allenic cycle, which is also predicted at much higher levels of theory.[30] The OK- and ONa-substituted molecules form nonplanar crowns of the D3 symmetry. The hypothetical Li and BeH, as well as OLi-, SLi-, OBe-, and SBe-substituted molecules, are starlike D6 structures with metal atoms placed at the CC bond intersection straight lines. Most molecules with bulky substituents, like ONO2, OPh, OSiH3, and SMe, have the most stable forms with the groups alternatingly directed up and down the ring and thus exhibit the S6 symmetry. The other molecules also adopt similar conformations and have obvious features like hydrogen bonding (COOH, OH, CHO, etc.).

Finally, it must be stressed that quite a large group of molecules gathered in Table have been known for a long time. Mellitic acid (hexacarboxybenzene) and mellitates have been known since the first half of the 19th century.[31,32] The hexa iodo-, bromo-, chloro-, hexamethyl-, and hexahydroxy-benzenes (and hexahydroxylates) have also been known since the 19th century.[33−40] Hexafluorobenzene was only first synthesized in 1956,[41] but now, it is widely used as a valuable solvent and standard in 19F NMR. Hexaaminobenzene was first synthesized in 1929, and its crystalline state was already characterized in 1931.[42,43] Hexaphenylbenzene has been known since 1930[44] and hexa(trifluoromethyl)benzene, since 1960.[45] Hexacyanobenzene was first synthesized and patented in 1963.[46,47] A synthesis of explosive and propellant hexanitrobenzene was published in 1979,[48] hexacetylenobenzene (hexaethynyl-benzene) in 1986,[49] and the first synthesis of benzenehexathiol in 1989,[50] although the hexathiolates were known earlier.[51,52]

Inspection of Table reveals that OCN, OC≡CH, OMgH, and OH hexasubstituted benzenes exhibit HOMA = 1.000. Thus, they are geometrically isoaromatic with benzene. However, they are not magnetically isoaromatic as their NICS(1)ZZ indices are basically different: −22.2, −22.0, −19.7, −20.2, and −29.2, respectively. Analogously, there are many other pairs of molecules which are geometrically, but not magnetically, isoaromatic. Among these are, for example, SH and CN (HOMA = 0.944; NICS(1)ZZ = −16.1 and −20.5 ppm, respectively) or CF3 and NC (HOMA = 0.951; NICS(1)ZZ = −22.8 and −18.7 ppm, respectively). Also, pairs of magnetically, but not geometrically, isoaromatic molecules can be found: OMe and OCHO (NICS(1)ZZ = −22.0 ppm and HOMA = 0.988 and 0.994) and O– and OF (NICS(1)ZZ = −17.3 ppm and HOMA = −2.106 and 0.854). We think that the differences in magnetic and geometric isoaromaticities may be interesting in the design of new materials.

The search for geometrically ″overaromatic″ molecules stems from the observation that the NICS(1) and NICS(0) of some very common molecules, with symmetry close to D6, indicate that their magnetic aromaticity is greater than that of benzene, whereas the geometric aromaticity HOMA index suggests the opposite. Indeed, the NICS(1) values of as many as 12 hexahomosubstituted molecules listed in Table , namely, OCaH, NH3+, OMe, OMgH, OLi, OCaF, OC≡CH, F, OCN, CHO, C6, and BeH(☆), are lower than −11.305 ppm estimated for benzene at the B3LYP/6-31G** level. Simultaneously, their HOMA values vary from −0.2 to 1.0 (C6 vs OCN, OC≡CH, and OMgH, respectively, Table ). If a similar comparison is made for the NICS(0) index,[53] more than 30 hexahomosubstituted benzenes have values lower than −9.846 ppm determined for benzene. There may be several reasons for errors and discrepancies in NICS values: insufficient basis set, improper computational method (DFT functional), dependence of magnetic and tensor properties on the position of the origin of coordinates, and so forth. However, we first calculated the NICS(1)ZZ index, which is nearly independent of the ring σ-orbitals.[19] This index is based solely on the ZZ entry of the shielding tensor corresponding to the virtual (dummy) atom placed 1 Å above the ring center. It seems to be a fair single-point characteristic of the magnetic aromaticity.[54]

Analysis of the NICS(1)ZZ values yields a different picture of aromaticity than that based on a comparison of the NICS(1) or NICS(0) values (Table ). This is in agreement with the findings for monosubstituted benzenes studied by Krygowski et al.[55] and Feixas et al.[56] All of the considered molecules, except the unusual ones (Scheme ), exhibit NICS(1)ZZ values greater than that of benzene. This means that all of these rings are magnetically less aromatic than that of benzene. Thus, the crucial assumption in HOMA definition that benzene is the most geometrically aromatic molecule is strongly justified. Furthermore, benzene could also be assumed to be the most magnetically aromatic species, if the odd BeH(☆) molecule could be ignored.

So, let us now focus on this hypothetical and exceptional BeH(☆) molecule. Observe that two types of hexaberylhydride and hexalithium benzenes are listed in Table . The two types are planar and have D6 symmetry, but in the (☆) molecules, the BeH and Li moieties are attached through the middle of the C=C bonds, while in the (*) ones, through the C atoms (Scheme ). Although, the (☆) structures seem to be unusual, all their fundamental frequencies are real at the applied level of theory. In contrast, the (*) forms have structures typical for hexasubstituted benzenes, but some of their frequencies are imaginary indicating that they are not minima at the potential energy surfaces. In fact, the (☆) BeH and Li structures have energies much lower than the (*) ones, be they calculated using the 6-31G** or the much larger aug-cc-pVTZ basis set. The (*) Li and BeH isomers do not exhibit exceptional aromaticity features (Table ) and would not be worth mentioning if the NICS(1)ZZ of the BeH(☆) structure were not extraordinary. In conclusion, we speculate that the unusual BeH(☆) structure, which in normal circumstances would violently decompose, could be reconsidered and reevaluated in a future material chemistry search for a sophisticated superaromatic substance existing in some atypical environment. Thus, (☆) BeH and Li structures cannot be recommended as the aromaticity standard and the benzene molecule remains the single candidate for such a standard.

Does the Kekulé Ring Moiety Exist?

In the mid-1860s, Kekulé proposed the cyclic structural formula of benzene with alternating single and double bonds between the C atoms.[57,58] The Kekulé model improperly suggested the D3 symmetry of benzene, but at the end of the 19th century, Thiele introduced the bond resonance concept which removed this misconception.[59] Nevertheless, the Kekulé structural formula has been and still is ubiquitous in chemical literature as it is handy and everyone is aware of its conventionality. Besides, the Kekulé benzene plays an important role in the HOMA aromaticity index, which is defined to be exactly 0 for this hypothetical form.[1−3] Yet, in general, HOMA is a multivalent function of n variables (n = 6 for the six-membered ring) and can be zero for (potentially) an infinite number of sets of six CC distances. Note that HOMA can be equal to 0 even for the D6 symmetry moiety if the six bonds equally differ from Ropt by ±α–1/2 (eq ).

In the original paper, discussing the parameters of the HOMA index,[60] the values of Ropt = 1.388 Å and α = 257.7 Å–2 were justified by the experimental data. Now, in the vast majority of articles, aromaticity is discussed based on computational data for which Ropt and α should be readjusted. The Ropt value is usually assumed to be the CC distance in benzene optimized at a given level of theory. However, determining α and R(C–C) and R(C=C) in the cyclohexatriene moiety needs more assumptions. Krygowski’s assumptions[60] lead to the following resultsTo find rings satisfying eq , a series of triply fused hexasubstituted benzenes possessing the C3 axis (TFB molecules, Scheme ) were calculated using α = 257.7 Å–2 and the B3LYP/6-31G** calculated Ropt,6-31G** = 1.39632 Å. To facilitate the discussion, the systematic TFB names are simplified to the names of the triply fused moiety. Thus, here, triphenylene and trinaphthylene are called benzene and naphthalene derivative, respectively.

Scheme 2

Sketch of the Considered Triply Fused Hexasubstituted Benzenes (TFB)

FS denotes an arbitrary fused system.

Note at the beginning that low TFB symmetries like C3 or C3 appear when the fused systems are bulky and interact with each other. As a result, the whole molecule can be distorted or even chiral. This does not influence the HOMA index, which can be calculated straightforwardly, but the NICS(1) group of indices split into two values characterizing two faces of the ring rather than the ring itself.[21]

Observe also that despite the C3 symmetry of the TFB moiety, the HOMA index can change continuously from nearly 1.0 to less than −1.0 (Tables and S1). For example, the central ring in tris(cyclopenteno) benzene has HOMA = 0.987 and is aromatic while that in tris(cyclobutadieno) benzene has HOMA = −1.251 and is antiaromatic.

Table 2

Parameters of the Central Ring in TFB Represented by Molecules Fused through the Vertical Bond on the Left-Hand Side of the Cells Are in the 2nd and 9th Columnsa

The systems calculated at the B3LYP/6-31G** level and for HOMA, EN, and GEO, α = 257.7 Å–2 and Ropt,6–31G** = 1.39632 Å were assumed. No 1 = No 85 in Table S1.

In this context, we scan the chemical realm far beyond the Mills–Nixon hypothesis,[61] approximately saying that a small ring fused with benzene promotes the single- and double-bond alternation. A closer inspection into the ring condensation reveals that such a fusion influences an interplay between the single- and double-well potentials in benzene in many different ways.[62−64] Yet, here we focus on the C3 symmetry fused systems yielding HOMA of the central ring close to 0. In terms of the HOMA index, this denotes the perfect bond alternation. However, a more detailed look into the systems with HOMA ≈ 0 reveals unobvious peculiarities.

Let “f” and “b” denote the “ring fusion” and the “bay” CC bonds in the central benzene ring of TFB, respectively (Scheme ). Tables and S1 indicate that probably none of the relations between Rf, Rb, and Ropt listed in Table can be ruled out a priori.

Table 3

Different Types of TFB Molecules. Rf and Rb Denote the “Ring Fusion” and the “Bay” CC Bond Lengths in the Central Benzene Ring

TFB type	subtype	relation	θfb = Rf – Rb
TFB(i)	(a)	Rf > Ropt > Rb	+	–
	(b)	Rf < Ropt < Rb	–
TFB(ii)	(a)	Rf > Rb > Ropt	+	+
	(b)	Rb > Rf > Ropt	–
TFB(iii)	(a)	Ropt > Rf > Rb	+	+
	(b)	Ropt >Rb > Rf	–

Thus, a TFB derivative can exhibit: D3 symmetry, alternating bond lengths, and HOMA of the central ring close to 0, but it can simultaneously belong to the TFB(ii) or TFB(iii) type (Table ). However, based on Krygowski’s assumptions (eq ), these types are forbidden for the model cyclohexatriene moiety, which must belong to the TFB(i) type and must satisfy that R(C–C) > Ropt>R(C=C) (Figure a).[60] Only in such a case will the τfb ratio be negative (Table ), and for the perfect cyclohexatriene match, it will be either −2 or −1/2 (eq ).

Figure 1

(a) Plot of Rb (“bay”) versus Rf (“ring fusion”) CC bond lengths in TFB molecules (represented by molecules fused through the vertical bond on the left-hand side of the structural formulas) indicating sections of allowed and forbidden models of cyclohexatriene structure. (b) Plot of the difference of EN and GEO structural co-indices of HOMA versus HOMA demonstrating that abscissa = −0.8 indicates structures corresponding to proper relations between Rb, Rf, and Ropt.

Surprisingly, stable and relatively unstrained TFB candidates for the cyclohexatriene moiety, in which triply condensed moieties are fused through the six-membered ring, belong to the TFB(ii) group. Indeed, the HOMA index of the central rings in: benzothiophene (1), phenanthrene (2), hydroxyindoles (3, 27), benzofuran (5), fluoroindoles (10, 27), nitroindoles (12, 38), methylindole (16), benzoimidazole (18), N-methylindole (20), indole (21), diborinine (23), diaminonaphthalene (31), dihydroxynaphthalene (32), di(dimetyloamino)naphthalene (34), difluoronaphthalene (37), dicyanonaphthalene (39), dichloronaphthalene (4), dimethylnaphthalene (41), phthalazine (42), naphthylborane (43), and naphthalene (44) TFB derivatives ranges from 0.1 to −0.1 (Table ).

Also unexpectedly, the TFB(i) molecules exhibiting HOMA ≈ 0 are fused to the central ring through the pentalene system (Table ). The nonalternant antiaromatic pentalenes with their 4n π electrons are reactive even at low temperatures.[65−67] Even if they were interesting for material design, they seem inappropriate as model reference molecules.

It is worth recalling that HOMA = 1 – EN – GEO[5,6] and that the EN and GEO indices can also be helpful in the search for the best TFB approximation of cyclohexatriene (Tables and S2). Note that EN and GEO factors, named after words “energetic” and “geometric”, respectively, are purely structural indices based solely on the CC distanceswhere Rav is the mean CC bond length and the other symbols are defined in eq .

Equations and 3 indicate that for cyclohexatriene, EN = 0.1 and GEO = 0.9, which provide an additional criterion for the search. However, alike Rf and Rb, the relationships between EN, GEO, and the molecular structure are implicit. This is because they are functions of the sum of the molecular distances (the average) and or the sum of their squares. Yet, juxtaposition of Figure a with the plot of EN-GEO against HOMA (Figure b) shows that structures with EN-GEO ≈ −0.8 correspond to allowed cyclohexatriene models (Rf > Ropt > Rb) while those close to EN-GEO ≈ 0.8 correspond to forbidden structures (Rf, Rb > Ropt). It is striking that desired TFB(i) molecules with negative τfb ratio have GEO, EN, and EN/GEO ratio close to that in cyclohexatriene (Table ). Thus, analysis of EN and GEO values can serve as alternative searching parameters, however, slightly less direct than the Rf, Rb, and Ropt distances.

This part of the study was done based on trial-and-error, as well as analogy, search methods. A more systematic and exhaustive investigation would provide more molecules satisfying the necessary condition for the best cyclohexatriene model. However, we do not have much hope that they will be simple and convenient.

This section demonstrates that finding a good model for the hypothetical cyclohexatriene moiety for which HOMA = 0 is difficult. Moreover, even if a candidate molecule has the HOMA index close enough to null, it often exhibits Rf and Rb bond lengths either both greater or both smaller than the optimal CC distance in the benzene molecule. In terms of EN and GEO, it has EN-GEO ≈ 0.8 (Table ). On the other hand, systems that satisfy the condition that Ropt distance is between Rf and Rb ones (or EN-GEO ≈ −0.8) are likely to be unstable and thus inconvenient as a reference molecule.

Thus, if we are interested in the HOMA definition rather than in a material design, it may be practical to refer in the HOMA definition not to the mysterious cyclohexatriene molecule, but to an existing one. The HOMA value of such a molecule could be determined with the same computational or experimental method as used for benzene. The α coefficient in eq could be determined based on the two standard molecules. We suggest that the cyclohexane chair conformer of the D3 symmetry,[68] easy to calculate at any theoretical level used to calculate HOMA, can be such a model molecule. The HOMA of the cyclohexane chair conformer is equal to −4.08706 at the B3LYP/6-31G** level (Ropt,6–31G** = 1.39632 Å and α = 257.7 Å–2). Assuming that this value is exactly −4.000 implies α = 253.29 Å–2 and HOMA of cyclohexatriene still equal to 0.000. Thus, it is pragmatic to

set HOMA equal to 1.000 for the CC distance in the benzene molecule of D6 symmetry,

set HOMA equal to −4.000 for the CC distance in the chair cyclohexane of D3 symmetry,

determine the α parameter based on (1), (2), and eq for the chosen level of calculations.

The above allows for conserving eq , but the α and Ropt parameters must be reestablished at each level of theory applied using the HOMA values of the chair cyclohexane of D3 symmetry equal to −4.000. Such an approach allows for obtaining congruent HOMA values at any computational or experimental level applied.

Can a Geometrical Aromaticity Index Be Defined Better Than HOMA?

In our earlier studies, we showed the HOMA index to be more than just a geometrical aromaticity index.[7,8] Indeed, the primary connection of HOMA with aromaticity is in choosing aromatic benzene as a reference for which HOMA equals 1. A secondary connection is through the normalization parameter making HOMA equal to 0 for the hypothetical, totally nonaromatic, Kekuléne form of benzene. However, using the HOMA index, we obtained valuable characteristics for both unsaturated and saturated structures, were they cyclic or acyclic.[7,8] Furthermore, we have shown the HOMA index to be a topological index in a kind of modified graph theory.[9] So what makes the HOMA definition universal enough to suit for the construction of different descriptors? This is because the HOMA index can be interpreted as a linear function of well-defined mathematical objects such as: (a) the square of a distance in an abstract n-dimensional space or (b) a variance, i.e., the second central moment.

(a) Consider an abstract n-dimensional molecular space, MS, in which a molecule is represented by a point X = (x1, x2, ... x), where the subsequent coordinates correspond to the selected CC bond lengths in this molecule. A comparison between two molecules, X and Y, can be expressed as a distance, d(X, Y) in the MS. In an n-dimensional Euclidean space, it is most natural to use the distance given by the square root of a sum of squares of the Cartesian coordinate differences: .

Such a comparison also makes sense if we are interested only in a difference between restricted fragments of the molecules, such as rings, to evaluate their aromaticity. Instead of discussing the details of the MS concept, here it is enough to state thatwhere d(·) is a distance in MS, X is a ring in a given molecule and Y = (y1, y2, ... y) stands for benzene in which all coordinates are equal to each other: y1 = y2 = ··· = y = y = d(CC), while const corresponds to α/n from eq . Additionally, if we assume that there exists an abstract reference molecule Y for which y1 = y2 = ··· = y = y = d(CCBenzene) but m ≠ 6, then we can calculate the HOMA index for any molecule composed of an arbitrary number of C atoms.

(b) In statistics, a kth central moment of probability function of a random variable x, μk, i.e., the expected deviation of kth degree from the expectation value, is defined as follows: μ⟨(x – ⟨x⟩)⟩ = 1/n∑(x–x̅), where the angle brackets ⟨•⟩ and x̅ denote the arithmetic mean. An ordinary kth moment, m=⟨x⟩ = 1/n∑(x), is not taken about the mean and .

Thus, if x is identified with the ith length R in a molecule (a ring) and with the optimal bond length Ropt, then HOMA is connected to the variance that is the second central moment for the multiple discrete random variables x1, x2, ... x, i.e., as followswhere α is as in eq . Note that eq allows n to be other than 6 because it can be understood as a number of different comparisons with the optimal bond length. So, the HOMA index can be calculated for any molecule composed of an arbitrary number of CC bonds.

Since HOMA can be expressed as a linear function of d2 or μ2, a better description of the geometrical aromaticity can be searched for among the extensions of the HOMA index, which is treated as an abstract distance[69] or as a statistical parameter.[70] Let us follow the second approach.

For the molecules considered in this study, supplemented by several simple molecules such as alkanes, alkenes, cycloalkanes, cycloalkenes, polyallenes, polyalkynes, their substituted derivatives, etc., let us plot μk central moments (k = 2, 3, 4 and μ2, μ3, and μ4 are variance, skewness, and kurtosis, respectively, taken about m1, i.e., the mean of CC distance in benzene) against the first central moment μ1 (Figure ). Note that the range of changes of μk decreases by 1 order from k to (k + 1). Observe also that there are always some scattered points over the expected regular polynomial relationship between μ1 and μk (Figure ).

Figure 2

Plot of the μ2, μ3, and μ4 central moments against the first central moment μ1, (a), (b), and (c), respectively, for the series of studied molecules at the B3LYP/6-31G** level. All values are multiplied 100 times.

Mark that a trivial geometrical index expressed as μ1 (x – x̅) is already bringing to light crucial similarities and differences among molecules or among their selected parts (Figure a). For molecules for which 100·μ2 varies from 0 to 5, 100·μ1 ranges from −20 to +20 (Figure a). However, to emphasize more subtle differences between the molecules, consider an index expressed as a sum of the first four central momentswhere and m1,ref is the mean in the reference molecule. Dividing μk by kth power of m1,ref guarantees the dimensionlessness of the Γ index.

Note that, for the reference molecule, Γ = 0. The majority of points in the plot of Γ against HOMA (Figure a) is placed on the branches of the inverse of a parabola as the largest Γ term is due to μ1. To better expose interrelations between HOMA and Γ, consider a geometrical auxiliary index, GAI, defined as followswhere the sgn(·) function assigns the value 1 to any positive number, 0 to zero, and −1 to any negative number.

Figure 3

(a) Relationship between Γ index constructed as a sum of the first four central moments μn (about m1,ref) and the HOMA index. (b) Relationship between the GAI and HOMA indices.

Remark that:

GAI = 0 for the reference benzene molecule,

GAI has neither a lower nor upper limit,

For the same HOMA, GAI can take values from positive, through zero, to negative,

Plot of GAI vs HOMA reveals different trends for different classes of molecules,

For HOMA ≈ 1.0, GAI ≈ 0.0 and GAI is likely to be useful as an auxiliary aromaticity index,

For HOMA ≪ 1.0, HOMA and GAI are general geometrical, rather than aromaticity, indices and are possibly useful for the characterization of different classes of molecules (Figure b).

The same can be said about the Γ index.

It is likely that the split of plots into series of regular tendencies reveals classes of differently strained molecules (Figure b). The upper boundary seems to be due to unstrained molecules (the black line going through the empty circles), while the lower one corresponds to linear polyunsaturated compounds such as polyallenes (the red line going through the red circles). Let us add that a slight nonlinear curvature of those trends results from eq , defining in fact a higher polynomial.

To illustrate that for molecules far from being aromatic (such as n-alkanes), the GAI and HOMA indices may be valuable geometrical descriptors, we correlated the boiling points of these molecules[71] with the two indices (Figure ). We used the extended all-transoid geometries of the most stable conformers. Yet, the number of alkane conformers increases exponentially with the alkane size, the intramolecular distances depend on conformation, and the distance-based geometrical indices vary. Therefore, to verify if the choice of the most stable conformers could provide fortunate but false trends, we estimated the studied n-alkanes conformational spaces at the molecular mechanics (MM) level. Next, all of the conformers were reoptimized with the B3LYP/6-31G** method. Then, the two indices were averaged using the Boltzmann population factors of the conformers[9]where INDEX stands here for HOMA or GAI, wi is the population factor of ith conformer calculated using either the total or the Gibbs free energy, m is the number of conformers, and the subscripts w and i stand for the population weighted and the current individual conformer, respectively. Note that eq defines the index which contains energetic factors and thus is not purely geometrical. Still, a purely geometrical weighting factor can be obtained by normalization against the sum of indices of the entire population.

Figure 4

Nonlinear correlations between boiling points of n-alkanes and (a) GAI and (b) HOMA indices obtained using geometry of the most stable conformers (blue) and indices averaged over all n-alkane conformers using weighting factors based on total energies (red) and Gibbs free energies (green).

It appeared that the correlations estimated for the most stable conformers have the same nonlinear and monotonic characters as those based on the Boltzmann population factors (Figure ). Hence, to check whether analogous relationships exist for arbitrary alkanes (up to nonane), we considered all of their constitutional isomers represented by MM-selected the most stable structures reoptimized at the B3LYP/6-31G** level. The plot of the so calculated HOMA index against the boiling points reveals more complex associations between the two variables (Figure ).

Figure 5

(a) Associations between the boiling points of n-alkanes and the HOMA index. (b) Nonlinear correlations between the boiling points of the most branched constitutional isomer of given alkane and its HOMA value. The HOMA values are based on structures of the most stable conformers. The same color in (a) and (b) denotes the same kind of alkane.

We observe that the HOMA value of the n-alkane is the greatest among the constitutional isomers (HOMA > −4.0, Figure a). The correlation observed in Figure can also be recognized on the right-hand side of Figure a. On the other hand, the most branched constitutional isomer has the smallest HOMA value (Figure a). The boiling points of the most branched constitutional isomers can also be nonlinearly correlated with HOMA (Figure b). Note that the boiling point of the most branched alkane is not always the lowest boiling point among the isomers (Table S2). Observe also that the boiling points of alkane constitutional isomers tend to increase with HOMA from the most branched to the most extended isomer. However, for lower alkanes, the correlation with HOMA seems to be strong (R = 0.95 for hexanes) while already for nonanes it is weak (R = 0.3) if it exists at all (Figure a). The analogous holds true for the GAI index, which for alkanes behaves quite similarly.

The connections between alkane branching and HOMA can be well explained. In the set of constitutional isomers, the n-alkanes are the most relaxed while the most branched constitutional isomers have the strongest intramolecular hindrances and internal repulsions. Therefore the average CC bond length in the former group is smaller than in the latter. In consequence, the difference between the CC distance in n-alkane and benzene is smaller than that for the most branched constitutional isomer. Hence, HOMA of the former is greater than HOMA of the latter.

However, the correlations and trends shown in Figures and 5 should be verified in a study devoted exclusively to this problem. First, the computational method should account for dispersion forces important for alkane intramolecular interactions and conformation. Second, the basis set used should be larger and should contain diffuse functions important to properly consider interactions between side chains. Third, the conformational landscape should be probed more carefully and should include methyl group rotations omitted here. Last but not least, several alkanes are chiral[70] and their number increases faster than the number of alkanes for which the optical isomerism is ignored.[72] Yet, if the number of asymmetric carbon atoms is greater than 1, the boiling points of different optical isomers may differ significantly. At the moment, the boiling points of chiral alkanes are determined only rarely[71] and the tabularized b.p. values correspond to racemic mixtures averaged over optical isomer populations. However, when simulating the conformational space, the chirality of alkanes must be taken into account and the derived parameters must be properly averaged. All of the above makes the study of flexible molecules quite a strenuous task that can hardly be accomplished in this study.

Summing up this section, there are at least two ways to obtain new geometrical indices: (a) by generalizing connections between the HOMA index and a distance in an abstract molecular space or (b) by generalizing the HOMA index connections to statistical central indices. Following the second way, we demonstrated that valuable indices, such as Γ or GAI, can be obtained. However, except for the fact that they have no upper and lower limits, the indices themselves are merely different from HOMA rather than better. Yet, juxtaposition of GAI and HOMA indices in common plot reveals additional molecular features such as molecular strain. However, a closer examination of this finding goes beyond the aim of this paper.

Conclusions

HOMA is one of the most popular indices for determining aromaticity and practically the only one which reflects the geometrical aspect of aromaticity. However, the HOMA index also performs well for any type of molecules, be they aromatic, only cyclic, or even linear. If so, it is important to better understand its behavior. To this aim, we have asked three questions which could provide a deeper insight into comprehension of this descriptor. This is also important because of the possible use of this simple index in molecular prescreening for further material design.

The questions are a bit provocative. The question “Does an “overaromatic″ carbon ring exist?” could not be answered based solely on the HOMA index, which inherently assumes that benzene is the most aromatic system. Therefore, the HOMA indices of hexahomosubstituted benzenes were juxtaposed with different magnetic aromaticity NICS indices, which have no such limitation. As a result, we find that although the less precise NICS(1) index would suggest that some of the considered molecules could be “overaromatic,” use of the NICSZZ(1) index almost unequivocally indicates that benzene deserves to be considered as the most aromatic molecule.

The answer to the question “Does the Kekulé ring moiety exist?” was based on analysis of triply fused hexasubstituted benzenes (TFB) with the C3 axis. Such molecules have two types of CC bonds in the central ring, “f” and “b”, corresponding to the “ring fusion” and the “bay” part of TFB, respectively. However, we demonstrated that a TFB derivative can exhibit high D3 symmetry, alternating bond lengths, and HOMA close to 0, but it can simultaneously have Rf, Rb, and Ropt both greater than or smaller than the Ropt of benzene. Yet, these types are forbidden to the model cyclohexatriene moiety, which must satisfy R(C–C) >Ropt > R(C=C). However, most of the molecules found which exhibited HOMA ≈ 0 and satisfied this condition contained the reactive pentalene system fused to the central ring. We also demonstrated that analysis of EN and GEO indices closely related to HOMA can provide the same conclusions. Hence, even if such THB molecules were interesting for material design, they present as inappropriate model reference molecules. We conclude that insofar as we are interested in the HOMA definition rather than material design, it could be practical to replace the mysterious cyclohexatriene molecule with an existing one. We suggest that the cyclohexane chair conformer of the D3 symmetry could be such a model molecule. The assumption that HOMA of the cyclohexane chair conformer is equal to exactly −4.000 implies α = 253.29 Å–2. This allows us to conserve the current definition, and HOMA of the Kekuléne cyclohexatriene is still equal to 0.000. But for each level of theory applied, it is necessary to reestablish α and Ropt parameters using HOMA of the chair cyclohexane equal to −4.000.

The third question is: “Can a geometrical aromaticity index be defined better then HOMA?” To find an answer to it, we have first demonstrated why the HOMA index is very successful. This is because it is connected to well-defined mathematical objects such as the distance in an abstract molecular space or the second central statistical moment. Based on the first four central statistical moments, we proposed a new geometrical auxiliary index (GAI) which has several properties similar to HOMA, but neither an upper or lower limit. As HOMA, it performs well both for aromatic and nonaromatic molecules as well as for cyclic and linear ones. Juxtaposition of the GAI and HOMA indices shows splitting of the plot into several regular tendencies, which ostensibly reveal classes of differently strained molecules. Moreover, GAI can probably serve as an auxiliary geometrical aromaticity index. Finally, we have illustrated that for molecules such as types of alkanes, both the GAI and HOMA indices may be valuable geometrical descriptors correlating with the alkanes boiling points. Nevertheless, for n-alkanes, the most branched alkanes one can find good correlations with boiling points, the general relation is complicated if it exists at all. However, further studies of this interesting feature go beyond the purposes of this study.