Polyradical Character of Triangular Non-Kekulé Structures, Zethrenes, p-Quinodimethane-Linked Bisphenalenyl, and the Clar Goblet in Comparison: An Extended Multireference Study

In this work, two different classes of polyaromatic hydrocarbon (PAH) systems have been investigated in order to characterize the amount of polyradical character and to localize the specific regions of chemical reactivity: (a) the non-Kekulé triangular structures phenalenyl, triangulene and a π-extended triangulene system with high-spin ground state and (b) PAHs based on zethrenes, p-quinodimethane-linked bisphenalenyl, and the Clar goblet containing varying polyradical character in their singlet ground state. The first class of structures already have open-shell character because of their high-spin ground state, which follows from the bonding pattern, whereas for the second class the open-shell character is generated either because of the competition between the closed-shell quinoid Kekulé and the open-shell singlet biradical resonance structures or the topology of the π-electron arrangement of the non-Kekulé form. High-level ab initio calculations based on multireference theory have been carried out to compute singlet–triplet splitting for the above-listed compounds and to provide insight into their chemical reactivity based on the polyradical character by means of unpaired densities. Unrestricted density functional theory and Hartree–Fock calculations have been performed for comparison also in order to obtain better insight into their applicability to these types of complicated radical systems.


INTRODUCTION
Over the past decade, graphene 1−3 has attracted considerable attention because of its wide range of applications, 4−9 e.g., as chemical sensors, organic semiconductors, energy storage devices, in spintronics and nonlinear optics. When the graphene sheet is cut into nanosized fragments, nanographenes containing several polycyclic aromatic hydrocarbon (PAH) units are formed. Therefore, the properties of PAHs are closely related to that of nanographene. These PAHs are composed of fused aromatic rings, a feature with almost unlimited possibilities that can lead to a rich diversity of compounds.
Many of the PAHs have a closed-shell electronic configuration in their ground state. However, there are types of PAHs which possess a high-spin, open-shell radical character in their ground state. 10,11 For example, phenalenyl 12 (1) (Chart 1) in its neutral ground state contains an odd number of carbon atoms with an odd number of π electrons, which makes it a radical. The extension of benzene rings in a triangular fashion can lead to several π-conjugated phenalenyl derivatives such as triangulene 13−16 (2), π-extended triangulene system 12,17 (3) (see Chart 1) and even larger systems.
These molecules are categorized as open-shell non-Kekuleṕ olynuclear benzenoid molecules where some electrons are unpaired due to the topology of the π electron arrangement. 18 The topological effect on the ground state of these molecules can be explained using Ovchinnikov's rule 18 which states that the ground-state spin quantum number, S, of a given PAH with six-membered rings can be expressed as S = (N* − N)/2 where N* and N are the number of starred and unstarred atoms of alternant π systems. As shown in Chart 1, in the case of phenalenyl (1) there are seven starred and six unstarred atoms. The open-shell character of this system can be understood qualitatively by considering no matter what resonance structure is chosen one starred atom will be left without a bonding partner. In the case of 2, two such empty valences are created, while there are three in the case of 3. Application of this rule leads to the aforementioned high-spin ground state. On the other hand, zethrenes 19,20 starting with heptazethrene 21,22 (4) (Chart 2), n-acenes longer than pentacene, 23−27 and some zigzag-edged graphene nanoribbons 27−31 (GNRs) have a singlet ground state but nevertheless significant biradical or even polyradical character. Among the different types of PAHs, zethrenes have attracted significant attention recently. 32−34 Zethrenes are z-shaped PAHs where two phenalenyl systems are head-to-head fused with or without a benzenoid spacer. The smallest member is the hexazethrene which contains a total of six condensed rings where the two phenalenyl rings are directly fused. The next one is heptazethrene (4) where an additional benzene ring is inserted in between the two phenalenyl units. If the terminal naphthalene units are replaced by anthracene, 34 depending on the position of the ring fusion, two different structural isomers with the same chemical compositions are formed: 1,2;9,10-dibenzoheptazethrene (5) and 5,6;13,14-dibenzoheptazethrene (6) (Chart 2).
As an alternative to zethrenes, Kubo and co-workers designed a p-quinodimethane-linked bisphenalenyl 35 (7) that contains p-quinodimethane and two phenalenyl rings (Chart 3). A characteristic feature of the zethrenes and structure 7 is the competition between a closed-shell quinoid Kekuléform and an open-shell biradical resonance form (Charts 2 and 3). 32 Because of the presence of these two resonance structures, two interesting questions arise: which resonance form is dominating in the ground state and, consequently, which isomer contains a larger biradical character? Clar's aromatic sextet rule 36−38 qualitatively predicts that for benzenoid PAHs the molecule with more aromatic sextet rings possesses more aromatic stabilization energy. Thus, if in different valence bond (VB) structures for a given PAH the biradical form contains more aromatic sextet rings than the closed-shell quinoid structure, then, as discussed in ref 32 and shown in Charts 2 and 3, its enhanced stability should be the source of a greater singlet biradical character. For example, in the case of structures 4 and 5 three aromatic sextets occur in the biradical form, whereas in the closed-shell Kekuléform only two aromatic sextets are present; for structure 6, a total of five sextets can be drawn (6c′) when the two radical centers are located at the terminal units. Thus, structures 4 and 5 are expected to exhibit significant biradical character, and the even larger number of sextets for structure 6 should help to stabilize the biradical form in this case even more. Chart 3 represents the resonance structures between Kekuléand the biradical forms for structure 7 where in the biradical form the central benzene ring obtains aromatic character. It is expected, therefore, that structure 7 also possesses a significant biradical character. Structure 8 in Chart 3 represents the non-Kekulébiradical form of the Clar goblet.
PAHs with radical character possess unique electronic, magnetic, and optical properties. Quantum chemical calculations play an important role for understanding the unique electronic characteristics of the PAHs. Among the methods available, density functional theory (DFT) would be the first choice because of its computational efficiency and good overall performance. However, due to the partly open-shell character of these compounds, an unrestricted approach has to be used which suffers from the problem of spin-contamination, 25,39 especially in the case of broken symmetry (BS) DFT for singlet biradical species or low-spin cases in general. Therefore, several other approaches such as spin-flip methods, 40,41 projected Hartree−Fock theory, 42 active-space variational two-electron reduced-density-matrix (2-RDM), 43,44 density matrix renormalization group (DMRG), 45,46 and coupled cluster with singles, doubles, and noniterative triples CCSD(T) 47,48 have been applied successfully to investigate the electronic structure of various PAHs. As a systematic and general alternative theory, multireference (MR) 49 methods have been used with great success as well. Among them, the multireference averaged quadratic coupled cluster (MR-AQCC) 50 approach represents a good option since it allows the inclusion of quasi-degenerate configurations in the reference wave function and of dynamic  electron correlation including size-extensivity contributions by  considering single and double excitations explicitly. Recently,  MR-AQCC calculations for the n-acenes 27,51,52 and periacenes 27,51,52 and several smaller challenging biradicals 53,54 have been performed successfully.
To understand the structural variation of the PAHs, their unique properties, and possible applications, it is important to investigate their electronic structure in more detail. Though Ovchinnikov's rule 18 and Clar's aromatic sextet rule 36−38 give a general idea about the multiplicity, stability, and electronic character of the ground state, it is nevertheless crucial to obtain a reliable quantitative description of these properties as well. Phenalenyl, its extended versions, and the combination to form zethrenes or p-quinodimethane-linked π-conjugated compounds such as the Clar goblet (Chart 1−Chart 3) constitute very interesting but also computationally challenging cases for describing the polyradical character of these compounds. Here, we use the aforementioned MR-AQCC method and the related MR configuration interaction with singles and doubles (MR-CISD) to discuss the properties of these types of PAHs in their ground state as well as in their lowest lying excited states. One major effort is dedicated to the clarification of which electronic configuration, either the closed-shell quinoid Kekuleó r the open-shell biradical form, better describe the ground state of structures 4−7. The related question of the multiplicity of the ground state and of singlet−triplet splitting will be addressed also by these high-level calculations. The complicated nature of the multireference wave function will be transformed to simple descriptors such as natural orbital (NO) occupation numbers and the unpaired density distribution which allow a concise chemical assessment of the polyradical nature of a compound and the location of chemically reactive unpaired density regions in a molecule. By means of comparison with our multireference ab initio results, DFT is also considered in the present context in order to obtain better insight into its applicability to these difficult questions concerning the correct description of biradical systems.

COMPUTATIONAL DETAILS
Multiconfiguration self-consistent field (MCSCF), mostly in the form of the complete active space self-consistent field (CASSCF) method, MR-CISD, and MR-AQCC calculations have been performed for the structures shown in Charts 1−3. For most of the calculations, two different sets of complete active spaces (CAS) have been chosen: (a) a CAS(7,7) and (b) a CAS (4,4). In active space (a), seven electrons have been distributed in seven orbitals and for active space (b) four electrons have been distributed over four orbitals. The choice of these active spaces is based on the NO occupation numbers computed from the unrestricted Hartree−Fock (UHF) wave function as suggested by Pulay et al. 55 The contribution of the single-and double-substituted configurations has been monitored and in the case of a weight larger than 1% extensions of the active space have been made in order to account for intruder states in the MR-AQCC calculations or to ensure degeneracies in the case of the triangular symmetry of structures 1−3, since for practical reasons only C 2v symmetry was employed. The active spaces for the MCSCF, MR-CISD, and MR-AQCC calculations were usually kept the same. Size extensivity corrections are computed for the MR-CISD approach by means of an extended Davidson correction, 49,56,57 denoted as +Q (MR-CISD+Q). More details on the construction of the MR wave functions and the orientation of the molecules in the Cartesian coordinate system can be found in the Supporting Information (SI).
Three different orbital freezing schemes have been used throughout the calculations: (i) core freezing where only the 1s core orbitals of all the carbon atoms are frozen at the MR-CISD level; (ii) σ-partial freezing (designated as (σ)+π-space) where 1s core orbitals of the carbon atoms (n) are frozen at the MCSCF level and depending on the degeneracy of the orbitals, additional n or (n − 1) σ orbitals, respectively, are frozen at the MR-CISD level; and (iii) π-space where all occupied and virtual σ orbitals were frozen by transforming the one-and twoelectron integrals from the atomic orbital (AO) basis into the basis of SCF orbitals keeping only the π orbitals active. The 6-31G, 6-31G*, and 6-311G(2d) 58,59 basis sets have been used throughout the calculations.
For the analysis of the radical character, we have used (a) NO occupations as computed from the AQCC density and (b) the unpaired density and the number of effectively unpaired electrons (N U ) 60−63 as originally introduced by Takatsuka et al. 60 as the distribution of "odd" electrons. This method was further developed by Staroverov and Davidson. 61 In order to emphasize the contributions from orbitals with occupation near one and suppressing contributions that are nearly doubly occupied or nearly unoccupied, the nonlinear formula of Head-Gordon 62,63 is used in which n i is the occupation of the ith NO and M is the total number of NOs. All the structures have been optimized using DFT with the B3LYP functional 64−66 and the 6-31G* basis set. Harmonic vibrational frequencies have been computed, and none of the structures except the 1 1 A g state of structure 5 shows imaginary frequencies in C 2h symmetry, which demonstrates that all other structures are minima. Structure 5 possesses two out-of-plane imaginary frequencies (b g and a u ) in C 2h symmetry because of the steric congestion between the hydrogen atom of the terminal benzene ring of the anthracene unit and the hydrogen atom of central benzene ring. The final optimized structure 5 has C 1 symmetry, and the corresponding frequency calculation confirmed that it is a minimum. It is important to note that the experimental structure 34 with the heptazethrene core also deviates from planarity with a torsional angle of 33.9°. The MR-AQCC calculation in C 1 symmetry would have been too costly. Since the difference between the planar C 2h and nonplanar C 1 structures in terms of singlet−triplet splitting and biradical character computed at DFT level is relatively small (which will discussed below) and also for comparison purpose with the other isomer (structure 6), we have decided to focus on the π conjugation and keep structure 5 planar by using C 2h symmetry.
A wave function stability analysis 67 of the Kohn−Sham determinant 68 has been performed for the optimized 1 1 A g state of structures 4−8. It was found that all the structures have a triplet instability and appropriate geometry reoptimizations have been performed in reduced planar C s symmetry. The real nature (except planar structure 5 has two imaginary frequencies) of the vibrational frequencies also confirmed that these structures corresponded to minima. The triplet states (1 3 B u for structures 4-6 and 1 3 B 2u for structures 7 and 8) have also been optimized separately. In the spin-unrestricted Kohn− Sham formalism used, especially the low-spin BS solution The Journal of Physical Chemistry A Article suffers from the problem of spin contamination. A high spin calculation with two unpaired electrons with parallel spins is applied to represent the triplet state, and the BS solution with antiparallel spins is used for the singlet state. To correct for spin contamination, the spin corrected formula as proposed by Yamaguchi 39,69 is employed. In this approach the vertical singlet−triplet gap at a given geometry "i" is given by the expression where E T i , E S i , and E BS i represent the energy of the triplet, singlet, and the BS solutions, respectively. ⟨S T i 2 ⟩ and ⟨S BS i 2 ⟩ are the expectation value of the square of the total spin operator for the triplet and the BS solutions.
For quantitative analysis of the open-shell character, along with N U , we have also computed the multiple diradical character 70,71 y i (i = 0, 1, 2) where 0 ≤ y i ≤ 1 and y i ≥ y i+1 from the spin-projected UHF (P-UHF) theory as where T i is the orbital overlap between the corresponding orbital pairs which can be expressed by where n HONO−i and n LUNO+i are the occupation numbers of the ith highest occupied NO (HONO) and the ith lowest unoccupied NO (LUNO) computed from UHF NOs. y 0 = 0 indicates pure closed-shell and y 0 = 1 indicates a pure diradical character. A perfect diradical has y 0 = 1 and y 1 = 0. Comparable y 0 and y 1 values indicate that in addition to the HONO/LUNO pair, non-HONO/LUNO pairs are also important. In the case of y 0 = 1 and a large y 1 value the diradical description is incomplete and a tetraradical character has to be considered. The MR calculations have been performed with the parallel version 72,73 of the COLUMBUS program system. 74−76 Population analysis of the unpaired densities has been carried out with the TheoDORE program. 77,78 For the DFT calculations along with the stability analysis and the UHF NO calculations, the TURBOMOLE program 79 has been used.

RESULTS AND DISCUSSION
3.1. Phenalenyl-Based Triangular Radicals. Figure 1 shows the MO occupation schemes for the ground state of the phenalenyl-based neutral radicals. The symmetry is given both in C 2v and D 3h (in parentheses) notations. The figure shows that the triangular PAHs exhibit a high-spin ground state and that the spin multiplicity increases with increasing molecular size. This is consistent with the non-zero value of (N* − N) as predicted by Ovchinnikov's rule 18 and also supports the experimental findings that the ground state of phenalenyl derivatives have a high-spin state. 12,14,16 In Table 1 a comparison of the energy differences between the ground and the degenerate first excited states of phenalenyl-based neutral radicals computed at MR-CISD and MR-CISD+Q levels is shown. Because of the severe occurrence of intruder states in the excited states, MR-AQCC calculations could not be performed in this case. Because of the use of the lower C 2v symmetry instead of the actual D 3h symmetry in the calculations, in some of the cases the degeneracy is slightly lifted. For this reason, we always considered the average energy between the two degenerate states. For phenalenyl (1), the degeneracy is well reproduced within ∼0.001 eV, and in most of the cases for the systems 2 and 3 it remains within ∼0.02 eV. The ground state of phenalenyl is 2 A 1 ″ (D 3h notation). The first excited doublet state is degenerate (E″ symmetry). The excitation energies amount to 3.165 eV at π-CAS(7,7) and 2.777 eV at π-MR-CISD levels, respectively, using the 6-311G(2d) basis and freezing all σ orbitals. Reducing the basis set to 6-31G* has only a minor effect (<0.02 eV). Inclusion of σ orbitals into the calculation increases the excitation energy at the MR-CISD level by ∼0.3 eV; the Davidson correction again decreases this value by 0.2 eV. At π-MR-CISD level, the occurrence of intruder states prohibited the use of the Davidson correction.
The next higher non-Kekule-type triangular PAH investigated is triangulene (2). Though it has an even number of carbon atoms, no Kekule-type formula can be given. 2,6,10-Tritert-butyltriangulene 14 is the first example of an experimentally a The symmetry is given as D 3h notation, D 3h (C 2v ) 1: . b For 1 and 3 a π-CAS(7,7) c and for 2 a π-CAS(4,7) c reference space are used, respectively. c Orbital occupation specifications are given in Table S1. d Contains intruder states.
The Journal of Physical Chemistry A Article synthesized genuine non-KekuléPAH. The experimentally observed linear dependence of the triplet signal intensity on 1/T for the 2,6,10-tri-tert-butyltriangulene combined with the observed g value of 2.0034 indicates that the high-spin triplet is the ground state. 14 The MR-CISD calculations collected in Table 1 for triangulene (2) confirm this finding. To characterize the manifold of lowest excited states for this system, the lowest triplet excitation energy and the triplet/singlet splitting are given in Table 1. The 1 3 E′ state is ∼2.7 eV higher than the 1 3 A 2 ′ ground state; the triplet/singlet splitting is only ∼0.8 eV.
Comparison of the excitation energies obtained with 6-31G* and 6-311G(2d) basis sets shows an agreement within a few hundredths of an eV. The same is true for the comparison with the 6-31G basis set (see Table S2). Freezing of the σ system also shows similar small energetic deviations.
Bearpark el al. 80 found for the trioxytriangulene trianion 16 at the CASSCF level the triplet state is more stable than the singlet state by 0.867 eV, which is in good agreement with our CASSCF results on triangulene (2). For the π-extended triangulene system (3) the ground state is a quartet, in agreement with Ovchinnikov's rule. At the π-MR-CISD level using the 6-311G(2d) basis, this state is 3.217 eV more stable than the first excited 4 E″ state, and the quartet/doublet splitting is 0.633 eV. The Davidson correction reduces these values by ∼0.2 eV. Again, the basis set effect is quite small. Analysis of the progression of the excitation energies along the increase of the triangular system from structures 1 to 3 shows relatively small effects. The excitation energy within the same spin multiplicity even increases from 2 to 3 by ∼0.3 eV, whereas the corresponding high spin/low spin excitation energies decreased by ∼0.15 eV.
The comparison of the energy differences for the open-shell triangular PAHs for the three different orbital schemes (core freezing, σ-partial freezing and π-space) using the 6-31G basis set is shown in Table S2. The differences in the results obtained with the polarized basis sets (Table 1) are mostly less than 0.1 eV. From these data, one can conclude that (i) π-space calculations are quite sufficient to describe the energy difference between the ground and the lowest excited states, and (ii) the basis set dependence on the addition of polarization functions is quite modest.
The NO occupations displayed in Figure 2 show that phenalenyl (1) has one (1a 1 ″), triangulene (2) has two (degenerate 4e″), and π-extended triangulene (3) has three (one in 2a 1 ″ and two in 6e″) singly occupied NO(s) (SONOs) in agreement with the MO occupation scheme given in Figure 1. It is important to note that the appearance of radical character with high spin multiplicities for the open-shell π-conjugated hydrocarbons arises because of the fused π-topology and not because of orbital degeneracy derived from the 3-fold symmetry. 12,70 Consequently, extending the π-conjugation according to the topology of π electrons can lead to an unlimited numbers of electron spins aligned parallel to each other in singly occupied MOs (SOMOs).
In Figure 3 the unpaired density is displayed for structures 1−3. It is delocalized over the entire molecule but resides mainly on the edges of the molecule. Moreover, for all three cases, the distribution at the edges is alternant and the unpaired density resides only on the starred atoms as indicated in Chart 1.
For phenalenyl, the central carbon atom does not carry any unpaired density. This depletion of unpaired density in the central region of the triangle has important consequences on the shape of stacked biradical dimers. For the phenalenyl dimer, it has been shown 81 that the dimer has a convex shape with the edge region being stronger bound due to the enhanced unpaired density in that region. Additionally, the possibilities for σand π-dimerization were discussed for a series of substituted phenalenyls based on a set of spectroscopic methods and on DFT calculations. 82 It is noteworthy that the total unpaired density values, N U presented in Figure 3, are larger as compared to the formal open-shell occupation of the high-spin state (1.33 vs 1 for phenalenyl, 2.50 vs 2.0 for triangulene, and 3.66 vs 3.0 for π-extended triangulene system). This excess unpaired density of these PAHs (1−3) is derived from nonnegligible contributions generated from NOs other than the SONO(s) (Figure 2   The Journal of Physical Chemistry A Article contribution to the singlet ground state. The difference in single to double bond lengths in the quinoid ring is between 0.06 and 0.08 Å for structure 4. A slightly lower range of 0.04−0.07 Å is found for both the C 2h -and C 1 -optimized structures of 1,2:9,10-dibenzoheptazethrene (5). This difference is significantly reduced for the alternative isomer, 5,6:13,14-dibenzoheptazethrene (6), to a range of 0.03−0.05 Å. This trend toward equidistant aromatic bond distances indicates the increasing importance of the biradical VB structure which will find its counterpart in the increase of the values for the total number of unpaired density (N U ) discussed below. The central benzene ring in structure 7 ( Figure S2) has almost aromatic character, indicating that this structure has a substantial biradical character. In the case of the structure of Clar's goblet ( Figure S2), all CC bonds are almost equivalent having a maximum difference of 0.046 Å with respect to the aromatic CC distance (1.39 Å). The exception is the central carbon−carbon bond which is strongly elongated (1.478 Å). This elongation divides the overall molecule into two separate moieties. The optimized bond distances a, b, and c (Chart 2) as shown in Figure S1  Singlet−triplet splitting is investigated next in order to continue the discussion of polyradical character of structures 4−8. In Table 2, computed vertical and adiabatic (in parentheses) singlet−triplet splitting energies ΔE(S−T) are listed and compared with available experimental data.
Vertical as well as the adiabatic ΔE(S−T) values for structures 4−8 are positive at all computational levels used, showing that these structures maintain a singlet ground state. However, there is a pronounced trend of ΔE(S−T) to decrease as one moves from structure 4 to 6. The relatively large adiabatic ΔE(S−T) for structure 5 (0.527 eV at π-MR-AQCC, Table 2) agrees well with the experimentally observed 34 sharp NMR peak and ESR silence. This means that the excited triplet state is located relatively high in energy. For structure 6, the adiabatic ΔE(S−T) is 0.324 eV (π-MR-AQCC), somewhat higher than that of the experimentally observed 32,34 value of 0.165 eV. Experimentally, 32 a broadening of the NMR spectrum at room temperature and a broad ESR signal were observed for structure 6 and the signal intensity is decreasing with decreasing temperature. This indicates the presence of a thermally excited  Table 2) showing that the singlet and the triplet states are almost degenerate. The evolution of the weights of the major configurations computed at π-MR-AQCC level (Table S4) goes in line with the just-described changes of the singlet−triplet splitting.
The singlet−triplet splitting has been computed also at the UDFT level and is given in Table 2 together with the expectation values of S 2 for the BS state. The spin contamination is quite significant for structures 5−8. Two kinds of ΔE(S−T) values have been computed at UDFT level: (i) using the spin-contaminated value and (ii) using the spin-corrected form (eq 3, ΔE proj ). The effect of the spin projection is to increase the UDFT singlet−triplet splitting, bringing it closer to the MR-AQCC values. For structure 5, ΔE proj of the nonplanar unrestricted C 1 structure is 0.568 eV, very close that of the value of 0.540 eV calculated for the planar C 2h structure.
A systematic comparison between singlet−triplet energies computed at several levels of orbital freezing schemes and basis sets including the unpolarized 6-31G basis (Table S5) finds in all cases a good agreement within a few tenths of an eV. These findings apply also to the computation of NO occupations (Table S6) and unpaired density (Table S7) discussed below. The observed relative insensitivity in terms of energies and character of the wave function facilitates the discussion of the polyradical character by means of multireference methods considerably. On the other hand, it is clear that this assessment has to be continuously re-evaluated when different kinds of π-conjugated systems are to be investigated.
NO occupations for the singlet state of structures 4−8 are presented in Figure 4. Several of the NO occupation numbers for the singlet state deviate strongly from the closed-shell limiting values of two/zero, indicating that they have substantial polyradical character.
Comparison of the HONO/LUNO occupation numbers (Table 3) indicates that deviation from the limiting values of two/zero is the smallest for structure 4 and the largest for structure 8. At the UHF level, these deviations from the closedshell reference values are considerably larger than the respective MR-AQCC results. In the former case, HONO/LUNO occupations are almost constant along the series 4−8, whereas at the MR-AQCC method there is a strong variation of the Table 2. Theoretical Vertical and Adiabatic (in Parentheses) Singlet−Triplet Splittings ΔE(S−T) (eV), Expectation Value of the Square of the Total Spin Operator for the BS Solution at the UDFT/B3LYP Level, and the Vertical Singlet−Triplet Gap ΔE proj Using the Spin-Corrected Formula for Heptazethrene (4), 1,2:9,10-Dibenzoheptazethrene (5), 5,6:13,14-Dibenzoheptazethrene (6), p-Quinodimethane-Linked Bisphenalenyl (7), and Clar Goblet (8)  The Journal of Physical Chemistry A Article occupation number, indicating a significant change in radical character. The picture of an almost uniform biradical character throughout the series 4−8 given by the UHF method is, however, not consistent with the graded evolution of the geometries and the singlet/triplet splitting discussed above. Even though most of the open-shell contributions computed at MR-AQCC level are coming from the HONO/LUNO occupation, for all structures, irrespective of singlet or triplet states, there are additional NOs, whose occupation numbers deviate significantly from the limiting value of two and zero ( Figure 4). This implies that in addition to the HONO/LUNO pair, other NOs also provide significant contributions to the radical character which cannot be neglected.
The densities of unpaired electrons for the singlet state are presented in Figure 5 and Figure 6 for structures 4−6 and 7−8, respectively. Unlike the situation found for the phenalenyl derivatives (1−3) where unpaired density is delocalized over the entire molecule, for structures 4, 5, and 6, the radical character is mostly distributed over a few positions. For structure 4, most of the unpaired density is located at C4/12 (see Chart 2 for numbering). For structure 5 the unpaired density is extended also to the atom pair C5/13. For both of these structures, the unpaired electron density within the benzene ring connecting the two phenalenyl segments is significant. In the case of structure 6, the unpaired density is strongly enhanced as compared to those of 4 and 5.
The main contributions are located equally at C4/12 and C7/15 indicating equal contributions from both the biradical resonance forms (6b′ and 6c′) as shown in Chart 2. However, the unpaired density situated on the other centers cannot be neglected. This indicates the existence of several additional VB structures in comparison to which are given in Charts 2 and 3.
The unpaired density for structure 7 ( Figure 6) shows a pattern which is more delocalized than the one indicated by the two mesomeric forms (7a′ and 7b′) given in Chart 3. Thus, in this case, the unpaired density can be better represented     The Journal of Physical Chemistry A Article by the resonance form 7c′, where the structure 7 can be considered as a combination of two phenalenyl systems linked by a benzene ring. In structure 8 (Figure 6), the unpaired character is mostly located at the zigzag edges with the largest contribution at their centers (position at C9/18, see Chart 3). It is also noted that for the zethrenes and the structures 5 and 6, the linking benzene ring seems to play a more important role (i.e., there is a significant amount of unpaired density located on this connecting ring relative to the total number of unpaired density) than in case of the vertical connections between subunits in structure 8. Table 4 compares the N U values computed at π-MR-AQCC level with the multiple diradical character indices, y i (i = 0,1,2..), obtained from P-UHF theory for structures 4−8. To obtain a common basis for comparison with the N U 's, the y i values were multiplied by a factor of 2. It is observed that for the singlet ground states of structures 4 and 5 the 2·y 0 values, which are computed from the HONO/LUNO UHF occupations, are almost twice of the N U values computed from the MR-AQCC HONO/LUNO occupations. Once the structures acquire more biradical character (for structures 6−8), the two values approach each other. This behavior is derived from the discrepancies in the NO occupation numbers computed with the two different methods ( Table 3). The y 0 values reported in ref 32 for the structures 4−6 are somewhat smaller as compared to our values. However,, the trend of increasing y 0 value as one moves from structure 4 to 6 is the same. For structure 5, the y 0 value as computed for C 1 optimized structure is 0.649, very close to that of the C 2h planar structure of 0.664. This indicates the similarity of the NO occupation numbers between the two structures. Comparing the N U values derived from different NO selections, it is noted that the total N U value is significantly larger than the one computed only from the HONO/LUNO part. These additional contributions come partly from the HONO-1/LUNO+1 set (Table 4), but also from the large number of NOs whose occupation numbers deviate from the 0/2 e occupations to a lesser extent. This is in contrast to the tetracyanoethylene anion dimer (TCNE 2 2− ) and neutral K 2 TCNE 2 system 83 where the effect of the non-HONO/LUNO pairs is practically negligible. In the present case the contribution of the non-HONO/LUNO pairs to the total N U value is almost 50% for the singlet biradicaloid structures of 4 and 5, but for the biradical structures 6-8, it decreases from 36% to 31%. Comparison of the N U values computed from the non-HONO/LUNO pairs for the singlet and the triplet states of structures 4-8 shows that they are almost identical. This indicates that the main difference is coming from the different occupations of the HONO and LUNO pair and the remaining contributions are quite the same. Table 4 shows that for all structures the triplet state maintains practically a constant N U value (from 2.589 e to 2.816 e), whereas for the singlet state of structures 4−8 a large change in the N U value (from 1.026 e to 2.880 e) is observed. As discussed just before, these differences come primarily from the HONO/LUNO pair (0.506 e to 1.974 e). For structures 4 and 5, this increase is only moderate from 1.026 e to 1.492 e, but from structures 5 to 6 it is relatively large (1.492 e to 2.241 e). This clearly indicates that strong variations in the polyradical character within the zethrenes can be achieved by means of relative modest changes in the π conjugation. For the singlet state of structure 7, the N U value is also very large (1.865 e) indicating significant singlet biradical character as well. Among the singlet state of all the structures, the Clar goblet (8) has the largest polyradical character.
Comparison of the distribution of unpaired densities between singlet and triplet states for structures 4−6 shows characteristic differences (Figures 5 and 7). These differences are naturally larger for the cases with smaller N U values in the singlet state (especially 4) since for the triplet state single occupation of the HONO/LUNO pair is enforced. For structure 4, the location of maximum density in both the singlet and the triplet states are same (C4/12, see Chart 2 for numbering). However, additionally, for the triplet state the unpaired density extends with significant populations on the C5/13 and C7/15 positions, respectively. Enhancement of similar atom position is also observed for the triplet state for structures 5 and 6. Even though the N U values between singlet and triplet start to come closer to each other, the weights on individual atoms still differ. For e.g. in the singlet state of structure 6, the unpaired density is equally distributed between C4/12 and C7/15 positions, respectively; for the triplet state the maximum of unpaired density is located at C7/15. This detailed insight into the unpaired density distribution should provide improved approaches to tune the singlet−triplet gap for these compounds. On the other hand, for the singlet and triplet states of both the structures 7 and 8 ( Figure 6 and Figure S3) the distribution of unpaired density is very similar in nature.
The Journal of Physical Chemistry A Article triangulene system (3) have been chosen. In the second case, a series of three zethrenes, heptazethrene (4), 1,2:9,10dibenzoheptazethrene (5), 5,6:13,14-dibenzoheptazethrene (6), and the p-quinodimethane-linked bisphenalenyl (7) have been investigated. Additionally, the non-KekuléClar goblet (8) has been studied. The motivation in choosing these two types of PAHs is that structures 1−3 already possess open-shell character because of their high-spin ground state, whereas for structures 4−7 the competition between a closed-shell quinoid Kekulévalence bond structure and an open-shell singlet biradical resonance form determines the actual electronic structure and the chemical reactivity. For structure 8, the topology of the π-electron arrangement of the non-Kekuléform is the characteristic feature. To get a reliable quantitative description of these interesting systems, high-level ab initio multireference approaches have been used. Unrestricted density functional theory and Hartree−Fock calculations have been performed for structures 4−8 also in order to assess their applicability to these molecular systems possessing a complicated electronic structure. The triangular structures 1−3 always have a nondegenerate high-spin state as ground state. The spin state increases with increasing molecular size as predicted by the Ovchinnikov's rule 18 and is also in agreement with ESR measurements of tri-tert-butyl-substituted phenalenyl, 12 triangulene, 14 and triangulene derivative. 16 The calculations also show that the lowest excited state is always degenerate. Although the unpaired density of the ground state of structures 1−3 is delocalized over the entire molecule, it mainly resides on one of the carbon sublattices, i.e., the starred atoms as defined above, and is for the most part located on the edges, independent of the size of the triangle. This localization of the chemical reactivity has important consequences on the lengths of the intermolecular CC bonds and the in general convex shape of stacked phenalenyl dimers as has been discussed in detail in ref 81. For the second class of systems (structures 4−8), the ground state is always singlet with a varying amount of biradicaloid (structures 4, 5), biradical (structure 7), or polyradical (structures 6 and 8) character. The triplet states of all of the structures have polyradical character. All indicators such as bond length alternation, singlet−triplet splitting, NO occupations, and unpaired densities clearly demonstrate that within the zethrene structure family the singlet state of structure 6 possesses a much larger polyradical character as compared to structures 4 and 5. Structure 7 also has open-shell singlet biradical character in its ground state. Interpretation of these results within the valence bond picture confirms that Clar's aromatic sextet rule can be successfully applied for the ground state of these types of systems, but for a concrete characterization of the chemical reactivity high level quantum chemical calculations are needed. Among structures 4−8, the Clar goblet (8) has the maximum polyradical character having nearly degenerate singlet and triplet states.
The low-spin broken symmetry state computed at the UDFT/B3LYP level is highly spin contaminated. Spinprojection increases the singlet−triplet gap and brings the DFT and MR-AQCC results into good agreement. However, NO occupations derived from UHF calculations in the spirit of the UNO-CAS method 55 show a strong overshooting of the deviations from closed-shell character for most of the singlet systems investigated and, as a consequence, also an overestimation of the polyradical character as measured by the y 0 and y 1 indices as compared to total numbers of unpaired electrons computed by the MR-AQCC method.
Analysis of our MR-AQCC results and also those of previous ones preformed on the singlet−triplet splitting in polyacenes 51 shows only a minor influence of basis set effects and of the amount of correlating σ orbitals. Though it is possible to perform large MR calculations by considering both σ and π electrons, this is an attempt to provide a guide for managing even larger systems by performing multireference correlation calculations for the π system only where such calculations including both σ and π electrons are too expensive.
In spite of the complicated structure of the multireference wave functions, the chemical analysis of the polyradical character is straightforward on the basis of the unpaired densities. Such an analysis is very helpful in locating the chemically reactive centers and indicating those regions on which to focus in order to stabilize the highly reactive polyradicals. It will also be possible to accurately assess the effects of different types of substituent attached to systems carrying polyradical character and provide pictorial information on concomitant changes in the chemical reactivity.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b12393.
Computational details, orbital occupation specification of all the structures, energy difference between the ground and the excited states with basis sets, natural orbital occupation and number of effectively unpaired electrons