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 follows  (1) where R i and R opt 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 Kekuleńe 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 1). 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

ACS Omega
Article (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 1) 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 D 6h 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 Among the molecules collected in Table 1, there are some exceptional cases (Scheme 1). The C 6 molecule converges to the D 3h allenic cycle, which is also predicted at much higher levels of theory. 30 The OK-and ONa-substituted molecules form nonplanar crowns of the D 3d symmetry. The hypothetical Li and BeH, as well as OLi-, SLi-, OBe-, and SBe-substituted molecules, are starlike D 6h structures with metal atoms placed at the CC bond intersection straight lines. Most molecules with bulky substituents, like ONO 2 , OPh, OSiH 3 , and SMe, have the most stable forms with the groups alternatingly directed up and down the ring and thus exhibit the S 6 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 1 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 hexahydroxybenzenes (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 19 F 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 1 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) Table 1 Notation of Substituents a a Virtual BeH and Li hexasubstituted D 6h benzenes. At the B3LYP/6-31G** level, the (*) types exhibit a few imaginary modes while the (☆) types show all real frequencies.

ACS Omega
Article ppm, respectively) or CF 3  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 D 6h , 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 1, namely, OCaH, NH 3 + , OMe, OMgH, OLi, OCaF, OCCH, F, OCN, CHO, C 6 , 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 (C 6 vs OCN, OCCH, and OMgH, respectively, Table 1). 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 1). 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 1), 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 1. The two types are planar and have D 6h 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 1). 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 augcc-pVTZ basis set. The (*) Li and BeH isomers do not exhibit exceptional aromaticity features (Table 1) 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 D 3h 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 Kekuleb enzene 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 D 6h symmetry moiety if the six bonds equally differ from R opt by ±α −1/2 (eq 1).
In the original paper, discussing the parameters of the HOMA index, 60 the values of R opt = 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 R opt and α should be readjusted. The R opt 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 results To find rings satisfying eq 2, a series of triply fused hexasubstituted benzenes possessing the C 3 axis (TFB molecules, Scheme 2) were calculated using α = 257.7 Å −2 and the B3LYP/6-31G** calculated R opt,6-31G** = 1.39632 Å.

Scheme 2. Sketch of the Considered Triply Fused
Hexasubstituted Benzenes (TFB) a a FS denotes an arbitrary fused system.

Article
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. Note at the beginning that low TFB symmetries like C 3 or C 3v 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 C 3 symmetry of the TFB moiety, the HOMA index can change continuously from nearly 1.0 to less than −1.0 (Tables 2 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.
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 C 3 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 1). Tables 2 and S1 indicate that probably none of the relations between R f , R b , and R opt listed in Table 3 can be ruled out a priori.

ACS Omega
Article can simultaneously belong to the TFB(ii) or TFB(iii) type (Table 3). However, based on Krygowski's assumptions (eq 2), these types are forbidden for the model cyclohexatriene moiety, which must belong to the TFB(i) type and must satisfy that R(C−C) > R opt > R(CC) (Figure 1a). 60 Only in such a case will the τ fb ratio be negative (Table 3), and for the perfect cyclohexatriene match, it will be either −2 or −1/2 (eq 2).
Also unexpectedly, the TFB(i) molecules exhibiting HOMA ≈ 0 are fused to the central ring through the pentalene system ( Table 2). 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 2  and S2). Note that EN and GEO factors, named after words "energetic" and "geometric", respectively, are purely structural indices based solely on the CC distances where R av is the mean CC bond length and the other symbols are defined in eq 1.
Equations 2 and 3 indicate that for cyclohexatriene, EN = 0.1 and GEO = 0.9, which provide an additional criterion for the search. However, alike R f and R b , 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 1a with the plot of EN-GEO against HOMA (Figure 1 b) shows that structures with EN-GEO ≈ −0.8 correspond to allowed cyclohexatriene models (R f > R opt > R b ) while those close to EN-GEO ≈ 0.8 correspond to forbidden structures (R f , R b > R opt ). 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 2). Thus, analysis of EN and GEO values can serve as alternative searching parameters, however, slightly less direct than the R f , R b , and R opt 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 R f and R b 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 2). On the other hand, systems that satisfy the condition that R opt distance is between Table 3. Different Types of TFB Molecules. R f and R b Denote the "Ring Fusion" and the "Bay" CC Bond Lengths in the Central Benzene Ring TFB type subtype relation

ACS Omega
Article R f and R b 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 1 could be determined based on the two standard molecules. We suggest that the cyclohexane chair conformer of the D 3d 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 (R opt,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 (1) set HOMA equal to 1.000 for the CC distance in the benzene molecule of D 6h symmetry, (2) set HOMA equal to −4.000 for the CC distance in the chair cyclohexane of D 3d symmetry, (3) determine the α parameter based on (1), (2), and eq 1 for the chosen level of calculations.
The above allows for conserving eq 1, but the α and R opt parameters must be reestablished at each level of theory applied using the HOMA values of the chair cyclohexane of D 3d 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, Kekuleńe 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 = (x 1 , x 2 , ... x n ), 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 that = − ·d X Y HOMA 1 const ( , ) 2 (4) where d(·) is a distance in MS, X is a ring in a given molecule and Y = (y 1 , y 2 , ... y n ) stands for benzene in which all coordinates are equal to each other: y 1 = y 2 = ··· = y n = y = d(CC), while const corresponds to α/n from eq 1. Additionally, if we assume that there exists an abstract reference molecule Y for which y 1 = y 2 = ··· = y m = y = d(CC Benzene ) 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: Thus, if x i is identified with the ith length R i in a molecule (a ring) and ̅ x with the optimal bond length R opt , then HOMA is connected to the variance that is the second central moment for the multiple discrete random variables x 1 , x 2 , ...
where α is as in eq 1. Note that eq 5 allows n to be other than 6 because it can be understood as a number of different

ACS Omega
Article 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 d 2 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 m 1 , i.e., the mean of CC distance in benzene) against the first central moment μ 1 (Figure 2). 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 2).
Mark that a trivial geometrical index expressed as μ 1 = (x i − x ̅ ) is already bringing to light crucial similarities and differences among molecules or among their selected parts (Figure 2a). For molecules for which 100·μ 2 varies from 0 to 5, 100·μ 1 ranges from −20 to +20 (Figure 2a). However, to emphasize more subtle differences between the molecules, consider an index expressed as a sum of the first four central moments Note that, for the reference molecule, Γ = 0. The majority of points in the plot of Γ against HOMA (Figure 3a) 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 follows = Γ ·Γ GAI sgn( ) 2 (7) where the sgn(·) function assigns the value 1 to any positive number, 0 to zero, and −1 to any negative number.
Remark that: 1. GAI = 0 for the reference benzene molecule, 2. GAI has neither a lower nor upper limit, 3. For the same HOMA, GAI can take values from positive, through zero, to negative,  It is likely that the split of plots into series of regular tendencies reveals classes of differently strained molecules (Figure 3b). 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 6, 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 4). 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, w i 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 8 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.
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 4). 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 5).
We observe that the HOMA value of the n-alkane is the greatest among the constitutional isomers (HOMA > −4.0, Figure 5a). The correlation observed in Figure 4 can also be recognized on the right-hand side of Figure 5a. On the other hand, the most branched constitutional isomer has the smallest HOMA value (Figure 5a). The boiling points of the most branched constitutional isomers can also be nonlinearly correlated with HOMA ( Figure 5b). 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  5a). 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 nalkanes 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 4  and 5 should be verified in a study devoted exclusively to this problem. First, the computational method should account for

ACS Omega
Article 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 NICS ZZ (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 C 3 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 D 3 symmetry, alternating bond lengths, and HOMA close to 0, but it can simultaneously have R f , R b , and R opt both greater than or smaller than the R opt of benzene. Yet, these types are forbidden to the model cyclohexatriene moiety, which must satisfy R(C−C) >R opt > 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 D 3d 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 Kekuleńe cyclohexatriene is still equal to 0.000. But for each level of theory applied, it is necessary to reestablish α and R opt 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.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsomega.9b02628. Table 2  The author declares no competing financial interest.

■ ACKNOWLEDGMENTS
Insightful and very constructive comments of anonymous referees to this paper, which allowed us to significantly improve presentation of the problems and clarity of the text, are very gratefully acknowledged. This work was supported by the National Science Centre in Poland Grant No. 2017/25/B/ ST5/02267. Calculations were performed at Sẃierk Computing Centre, National Centre for Nuclear Research, Sẃierk, Poland. A critical reading of the manuscript by Dr. Piotr F. J. Lipinśki from the Mossakowski Medical Research Centre, Polish Academy of Sciences, Warsaw, is gratefully acknowledged. The author thanks Rachel Muracka for her help in English language corrections.