Applications of H\"uckel-Su-Schrieer-Heeger method. I. Carbon-carbon bond lengths in polycyclic aromatic hydrocarbons

The equilibrium carbon-carbon bond lengths in $\pi$-electron hydrocarbons are very sensitive to the electronic ground-state characteristic. In the recent two papers by Stolarczyk and Krygowski (to appear soon) a simple quantum approach, the Augmented H\"uckel Molecular Orbital (AugHMO) model, is proposed for the qualitative, as well as quantitative, study of this phenomenon. The simplest realization of the AugHMO model is the H\"uckel-Su-Schrieffer-Heeger (HSSH) method, in which the resonance integral $\beta$ of the HMO model is a linear function the bond length. In the present paper the HSSH method is applied in a study of carbon-carbon bond lengths in a set of 34 selected policyclic aromatic hydrocarbons (PAHs). This is exactly the set of molecules analyzed by Rieger and M\"ullen (\textit{J. Phys. Org. Chem.} \textbf{2010}, \textit{23}, 315) in the context of their electronic-excitation spectra. These PAHs have been obtained by chemical synthesis, but in most cases no diffraction data (by X-rays or neutrons) of sufficient quality is available to provide us with their geometry. On the other hand, these PAHs are rather big (up to 96 carbon atoms), and \textit{ab initio} methods of quantum chemistry are too expensive for reliable geometry optimization. That makes the HSSH method a very attractive alternative. Our HSSH calculations uncover a modular architecture of certain classes of PAHs. For the studied molecules (and their fragments -- modules) we calculate the values of the aromaticity index HOMA.


Introduction
The variability of the (equilibrium) carbon-carbon (C-C) bond lengths in π-electron hydrocarbons and all-carbon molecules (fullerenes, nanotubes, and graphene) is a consequence of the coupling between the nuclear framework and the system of mobile electrons that occupy the molecular orbitals of the π-symmetry. This phenomenon poses a challenge to both experiment and theory.
On the experimental side one faces difficulties in obtaining good-quality monocrystals. Moreover, the usual diffraction techniques (by X-rays or neutrons) suffer from the limitations of the theoretical model which translates the diffraction data into the molecular geometry. This model, in which spherical atoms are subject to uncorrelated thermal motions, neglects the changes of the electronic density due to chemical bonding (X-ray diffraction), as well as the rigid-body motions of the whole molecule. The latter effects, affecting both the X-ray and neutron studies, are especially important for large and rigid structures of π-electron hydrocarbons.
On the theory side there are convergence problems of quantum-mechanical calculations with respect to the orbital basis sets and the electronic-correlation contributions. These problems become especially acute with the growing size of the studied molecules. As it has been already in the past, a carefully crafted semi-empirical quantum approach may come to the rescue.
In the two recent papers by Stolarczyk and Krygowski [1,2] a very simple Augmented Hückel Molecular Orbital (AugHMO) model has been designed for this purpose. It allows for a qualitative, as well as quantitative, study of the bond-length variation in general π-electron molecules (i.e., those which are planar or locally planar). The simplest realization of the AugHMO model, the the Hückel-Su-Schrieffer-Heeger (HSSH) method, was parametrized [2] for the π-electron hydrocarbons and all-carbon molecules. It was demonstrated [2] that the HSSH method is capable of a very good description of the C-C bond lengths in a variety of π-electron hydrocarbons and carbon systems, including fullerene C 60 , polyacetylene, and graphene.
In the present paper we employ the HSSH method to calculate the equilibrium C-C bond lengths for a set of the polycyclic aromatic hydrocarbons (PAHs). These are precisely the molecules whose spectroscopic properties were summarized and analysed in a review paper by Riegel and Müllen [3]. These PAHs have been obtained by chemical synthesis (many of them by Müllen and his coworkers), but in most cases no diffraction data of sufficient quality is available to provide us with their geometries. On the other hand, these PAHs are rather big (up to 96 C atoms), and ab initio methods of quantum chemistry are too expensive for a reliable geometry optimization.

HSSH calculations
The HSSH bond-length optimization procedure is a (simplified) analog of the molecular-geometry optimization techniques that employ the analytical gradients and Hessians of the total molecular energy [4]. In the HSSH method the hydrocarbon under study is fully characterized by its topology of C-C bonds (which is coded in the form of the topological matrix [1,2]).
At the start of the calculations all the C-C bonds are put equal to the value for benzene, R e ben = 1.397 Å. Then the usual Hückel calculations follow, yielding the molecular orbitals (MOs) and the πelectron bond orders p bond . This is the 0th iteraction, fully equivalent to the standard Hückel approach. In each subsequent step (iteration) the condition of the vanishing gradient of the total HSSH energy is enforced: this amounts to iterating the linear bond-order bond-length (BO-BL) relationship: where R e bond and p e bond , respectively, are the equilibrium bond length and the π-electron bond order for a given C-C bond. R o = 1.523 Å and x = 0.189 Å are two geometrical parameters of the AugHMO model of π-electron hydrocarbons [1,2]. In the HSSH method for hydrocarbons [2] there is also a third geometrical parameter y = 0.2756 Å which determines the slope of a linear function β(R) describing the dependence of the "resonance integral" β on the C-C bond length R. In order to improve the rate of convergence of the iterative procedure, the Hessian of the HSSH total energy is calculated and applied [1]. In our calculations for PAHs we assumed that the iterations stopped when the subsequent values of bond-length differences dropped below 5 · 10 −6 Å|. To this end five "direct iterations" using Eq. (1) followed by two to three iterations involving the Hessian were always sufficient.
Four points concerning the HSSH method [1,2] should be made clear. First: this method optimizes only the bond lengths, and it is assumed that the optimal valence angles (the C-C-C ones in the case of PAHs) are consistent with the optimized C-C bond lengths, and as close to 120 o as possible. By reducing the molecular geometry to bond lengths only seems justified for such molecules as PAHs and contribute to the computational effectiveness of the HSSH method. Second: while Eq. (1) looks like a relationship of local character, it should be remembered that the π-electron bond orders p bond are derived from the set of completely delocalized occupied π molecular orbitals [1]. And in the HMO model (as well as in the AugHMO one) these bond orders appear to be sensitive probes of molecular topology. Third: PAHs belong to the class of alternant π-electron hydrocarbons [2], and within the HMO (AugHMO) model the Coulson-Rushbrooke theorem [5] holds: (i) the π-electron orbital energies for the occupied and unoccupied states are placed symmetrically with respect to the value of the "Coulomb integral" α, (ii) the net π-electron charges [1] at the C atoms are equal to zero. Four: in the HSSH calculations for hydrocarbons we employ the Hückel energy units: where R e ben = 1.397 Å serves as the reference C-C bond length [2].Thus, the only empirical parameters of the HSSH model are the above-mentioned three geometrical parameters: R o , x, and y, plus R ben .
We used our HSSH optimized bond lengths to calculate the values of the aromaticity index HOMA of Krygowski and coworkers [6][7][8][9]. HOMA (Harmonic-Oscillator Measure of Aromaticity) is currently considered the most important indicator of the aromatic character of π-electron molecules (or their fragments), based solely on the values of bond lengths. For a π-electron hydrocarbon the definition of HOMA reads as [6,7] where the summation goes over all N carbon-carbon bonds in the π-electron hydrocarbon (or its fragment), α = 257.7 Å −2 is a normalization constant, while R opt = 1.388 Å represents the optimal value of the C-C bond length in the HOMA model. Although of energetic provenience [6], HOMA of Eq. (3) is dimensionless and fulfils condition 0 ≤ HOMA ≤ 1 (smaller HOMA indicates lower aromatic character). It was found [8,9] that definition (3) can be rewritten in a more revealing form: where is the average C-C bond length in the π-electron hydrocarbon (or its fragment). The components EN and GEO may be interpteted as some "dearomatization" contributions originating from the bondlength elongation and bond alternation, respectively [9]. Let us remark that the EN component is simply a quadratic function of R av , and that the EN and GEO components (and thus HOMA) are sensitive to the rounding errors in the values of R e bond .

Modular architecture of PAHs
The carbon skeletons of PAHs are finite jigsaw cuts from the honeycomb lattice of graphene. At the perimeter of a PAH there is a number of methine groups (>CH) corresponding to the tertiary carbon atoms (each with two C neighbors); the remaining C atoms are of the quaternary kind (with three C neighbors). The varied shapes of PAHs determine the corresponding topologies of the C-C bonds, which may be represented by the so-called Hückel graphs (see Ref. [10], p. 28). we shall use such graphs in the figures throughout the paper: the C atoms are graph vertices, and the C-C σ bonds are graph edges (the H atoms and the C-H bonds are absent). In the AugHMO model [1,2], the topologies of the C-C bonds directly translate into the bond-length pattern. Our HSSH calculations for PAHs reveal a considerable spread (ca. 0.1 Å) of the equilibrium C-C bond lengths: from 1.363 Å [bond a in pyrene (15), see Fig. 4], to 1.461 Å [bond m' in perylene (26), see Fig. 7]. These values are to be compared with the HSSH result for graphene [2] equal to 1.424 Å.
Our HSSH calculations indicate that some classes of PAHs (clarenes, rylenes, and a subclass of the K-region PAHs) are well represented by a modular architecture. The ground-state geometry of such a PAH can be assemblied by using certain standard molecular fragments (modules) of fixed geometries, connected by C-C bonds (linkers) of fixed lengths. The set of modules employed in the present study is depicted, in the form of Hückel graphs, in Fig. 1. These modules are assumed to have "natural" symmetries, independent of the particular molecular environment. More details will be provided in sections devoted to clarenes (Sec. 4), K-region PAHs (2) (Subsec. 5.2), and rylenes (Sec. 7). How the bond lengths corresponding to the modules and linkers were computed? Take the example of bond a in module A3: the length R e a was calculated by averaging over all the 28 symmetry-unrelated occurences of this type of bond in the molecules containig the A3-module (see Figs. 2, 3, and 5). As a measure of the spread of the actual values of this bond length we used the absolute value of maximal difference, corresponding to the above-mentioned set of 28 values of R e a .

Clarenes
Eric Clar proposed [11,12] that the all-benzenoid PAHs are special because of being ideal superpositions of six-electron units -the "Clar sextets." Thus, it seems to be quite proper to call this class of PAHs "clarenes." The family of clarenes, as taken from the paper by Riegel and Müllen [3], consists of 14 molecules, form convenience partitioned into two subsets shown in Figs. 2 and 3.  Every molecule of clarene contain 6n carbon atoms and can be partitioned into six-carbon-atom benzene-like fragments, hereafter called the A-modules. Each A-module carries six π electrons -this is the above-mentioned Clar sextet. Within a given molecule the A-modules (depicted in Fig. 1) are connected by C-C bonds, hereafter called the AA-linkers. The An-module (n = 2, 3, 4, 5 ,6) uses n linkers to connect with other modules; benzene may be considered a "honorary" A0 module, while the phenyl group (not appearing in PAHs) is the A1-module. Modules A4 and A4' correspond to two diffrent arrangements of four linkers. Moreover, linkers connecting the same pair of modules (eg., A3-A5) may correspond to different linking topologies: see, e.g., linker u in molecule 12, and linkers u' and u" in molecule 6.
Our HSSH calculations demonstrate that, with sufficient accuracy, each clarene can be additively assemblied from standard A-modules, connected via some standard AA-linkers. The results in Tabl. 1 and 2 are quite convincing: the An-modules (n = 2 -5) have well-defined geometries (with characteristic bond lengths), and the AA-linkers of a given type have nearly constant lengths. It is seen that within the A-modules the bond lengths are distictly shorter than those corresponding to the AA-linkers. Clearly, this is an indication that the Clar sextets are to some degree localized within the clarene molecule.

Module
Benz.
The An-modules (n = 2 -5) necessarily have to be located at the perimeter of a given clarene. The interior of the molecule, if sufficiently big, is filled with the A6-modules. Quite surprisingly, A6 cannot be considered a rigid building block: in a less symmetric surrounding its perfect 6mm (C 6v ) symmetry 1 becomes visibly perturbed. The most important cases of A6 deformations are presented in Tabl. 3. In molecules 11 and 20 (see Fig. 5) the molecular symmetry group 3m (C 3v ) induces some bond alternation in the central A6-module; however, this effect is absent in molecule 14 due to a different orientation of the molecular-symmetry elements with respect to the central A6 hexagon. A similar pattern of deformation is found in molecule 19, despite its lower symmetry corresponding to the m (C s ) group. It is apparent that bond alternation is the softest mode of deformation for the A6-module.  Table 3: A6-modules, HSSH equilibrium C-C bond lengths (in Å). In molecules 11, 19, and 20 a pattern of alternating bonds (a' a) 3 is seen. Full averaging over the remaining molecules containing A6-modules gives R e a = R e a = 1.417, ∆ = 0.003.
While it seems proven that the perimeter of a clarene assumes a definite geometric structure 1 We use here the symbols corresponding to the planar point groups corresponding to some arrangement od the An-modules (n = 2 -5), its interior is expected to approach the structure of graphene (however, see the discussion in the last Section). To this end the equilibrium bond lengths corresponding to bonds a and a' in the A6 module (see Fig. 1), and the equilibrium bond length corresponding to the A6-A6 linker (z", see molecule 14 in Fig. 3) should approach the HSSH value for graphene (1.424 Å). For the largest PAH in this study, molecule 14 (96 C atoms), we found R e a = R e a = 1.419 Å, R e z = 1.435 Å, quite far from the graphene limit. However, when one averages over three bond lengths adjacent to any C atom in any A6 module, the result is very close to 1.424 Å. Thus, even in relatively small clarenes their interiors "prepare" for becoming graphenelike. It may be said with some exaggeration that a big clarene molecule has a "hard shell" made of some An-modules (n = 2 -5), protecting an interior consisting of "softer" A6-modules, which, with the growing size of the molecule, loose their identity and melt into the homogenous honeycomb lattice of graphene.
The inspection of the HOMA values and its components in Tabl. 4 indicates that the An-modules retain, to some degree, the aromaticity of the benzene molecule. It is seen that it is EN (and thus R av ) which is mostly responsible for diminishing of HOMA with the increase of n (the number of linkers). The result for graphene suggests its low aromaticity, despite the fact that graphene is the most stable of all carbon and π-electron hydrocarbon molecules (as the molecular enthalpies of formation attest).  The molecular values of HOMA and its components collected in Tabl. 5 show that the aromatic character of clarenes is lowering with their size. Both the EN and GEO components contribute to this effect. As clarenes grow larger, one finds that (i) the value of R av increases and thus departs from the HOMA value of R opt , and (ii) the HOMO-LUMO gap (as calculated with the HSSH model) diminishes. Even if big clarenes are still far from the graphene limit, one may speculate that the metallic character and aromaticity are in conflict. It should be also remembered that in a PAH all the quaternary C atoms are completely inert [13], and thus all the chemistry takes place at the perimeter of the molecule, involving the methine groups. Therefore, the interior of PAHs (and clarenes in particular) may be of little (or no at all) relevance to their chemical reactivity.

K-region PAHs
Riegel and Müllen [3] coined the name "K-region PAH" to denote the PAH which can be derived from a clarene having one or more "bays" within its perimeter. The "bay" consist of four C atoms arranged as in the s-cis butadiene molecule (at the "mouth" of the bay there are two methine groups, causing some steric tension). Examples of such bays are seen in clarenes 1 -12 and 14. The K-region emerges when the mouth of the bay is closed with a two-carbon insertion of the formula -CH-CH-; we shall call this insertion the "f-handle", alluding to the bond f in phenanthrene (21), see Fig. 6. Obviously, the new molecule does not fit the clarene paradigm, and the K-region (corresponding to the f-handle) introduces some olefinic character to the molecule (for a detailed discussion on the chemical and spectral properties of the K-region PAHs, see Ref. [3]) The K-region PAHs of Ref. [3] may be conveniently divided into two groups: one containing pyrene and coronene (and ovalene, added here), and second composed of molecules which may be derived from "superbenzene" (9), see Fig. 3.

K-region PAHs (1): pyrene, coronene, and ovalene
Pyrene, coronene, and ovalene, depicted in Fig. 4, are well-known PAHs. Their C-C bond lengths are listed in Tabl. 6: a closer inspection of these data indicates that no obvious transferable fragments (modules) can be identified. This is a striking contrast with the clarenes of the previous Section. One can identify the "olefinic" f-handles as the a-bonds in all three molecules, and the f-bond in ovalene.   Table 6: Pyrene, coronene, and ovalene. HSSH equilibrium C-C bond lengths (in Å).
The molecular values of HOMA and its components are collected in Tabl. 7. Although the general trend is similar to that for clarenes, Tabl. 5, this time also the GEO contribution is important. Let us note that HOMA for coronene is much lower than for benzene, and not very different from that of graphene, see Tabl. 4.
In molecules (17 -20) one may identify the A3 and A6 modules corresponding to clarene 9, with practically "frozen" structures, see Tabl. 2 (also linkers r and v are unchanged). But a new module can be identified: this is the C4-module (see Fig. 1), involving fourteen C atoms (and fourteen π electrons), which resembles the phenanthrene molecule (9), the "C0-module". As seen in Tabl. 8, the (averaged) equilibrium C-C bond lengths of the C4-module are pretty well transferable, and so are the corresponding linkers (linker k being here a bit of exception).
The values of HOMA and its components for phenanthrene and the C4-module are collected in Tabl. 9. The molecular values of HOMA and its components presented in Tabl. 10 show that each insertion of the f-handle into clarene lowers its HOMA and thus diminishes its aromatic character.

Phenacenes
The first five members of the family of phenacenes are shown in Fig. 6. With the growing size they approach a limit of a π-electron polymer (not the 2D structure of graphene). The bond lengths of phenacenes are presented in Tabl. 11: here also we see no clue to identify transferable modules, although some regularities in bond lengths are clearly visible. Our HSSH equilibrium C-C bond lengths in phenacenes may be compared to those calculated by Firouzi and Zahedi [15] by means of the ab initio B3LYP/6-31G(d) approach. The agreement is satisfactory, with the exception of the "bridge bonds" (of the symbols d, h, and j), for which the HSSH values are smaller by ca. 0.02 Å (see also acenes, Sec. 8).
The molecular values of HOMA and its components for phenacenes are presented in Tabl. 12. It seems that these quantities may converge to some limiting values in the case of the polyphenacene polymer.

Rylenes
in Fig. 7 we present the first four members of the family of rylenes. As in the case of the family of phenacenes, with the growing size rylenes approach a limit of a π-electron polymer. However, in difference to phenacenes, rylenes display modular architecture, based on the ten C atom, ten πelectron B2 and B4 modules, see Fig. 1. They may be considered certain analogs of naphthalene (the "B0 module", see Fig. 8). The structures of the B2 and B4 modules are presented in Tabl. 13, together with the lengths of the corresponding BB-linkers.

Acenes
The first six members of the family of acenes are shown in Fig. 8. Acenes are similar to phenacenes of Sec. 6. They also, with the growing size, approach a limit of a π-electron polymer. The bond lengths of acenes are presented in Tabl. 16: here also we see no clue to identify transferable modules despite some regularities in bond lengths.  Table 16: Acenes. HSSH equilibrium C-C bond lengths (in Å).
We compared our HSSH equilibrium C-C bond lengths in acenes to those calculated by Firouzi and Zahedi [15] by means of the ab initio B3LYP/6-31G(d) approach. As in the case of the phenacenes, see Sec. 6, the agreement is satisfactory, with the exception of the "bridge bonds" (of the symbols d, g, and j), for which the HSSH values are smaller by ca. 0.02 Å.
The molecular values of HOMA and its components for acenes are presented in Tabl. 17. As seen, for long acenes HOMA quickly deteriorates, which is consistent with a high chemical reactivity of these molecules.

Concluding remarks
In a series of papers Tyutyulkov, Müllen, and their collaborators [16][17][18][19] posed a question: "is graphene an ultimate large hydrocarbon"? In order to find the answer, they analyzed various aspects: structural (topology), energetic (the energy spectra), and the influence of defects (including these with unpaired spins) and different edge structures. Their findings suggest [16] that the electronic-correlation effects should lead to a nonzero gap in the energy spectrum of PAHs, even in the limit of an infinite molecule (thus, strictly speaking, the metallic limit is not reached). However, the experimental energy gap in the clarene family (1 -14), as presented in Fig. 7 of Ref. [3], indicates that the band gap diminishes quickly with the size of clarene. Our HSSH study of the equilibrium C-C bond lengths for the PAHs considered in the review paper by Riegel and Müllen [3] is too limited to attempt an answer to the question of Refs. [16][17][18][19]. However, our elucidation of the modular architecture of some classes of PAHs (clarenes, rylenes, certain subclasses of the K-region PAHs) provides a suggestive indication that the properties of big PAHs may be understood on the basis of the calculations for smaller systems.