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 ) 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 ) and even larger systems.

Chart 1

Non-Kekulé Phenalenyl Based Neutral Radicals with Triangular Structures: Phenalenyl (1), Triangulene (2), and π-Extended Triangulene System (3)a

N* and N are the number of starred and unstarred atoms.

These molecules are categorized as open-shell non-Kekulé polynuclear 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 , 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 ), 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 ).

Chart 2

Structures of Heptazethrene (4), 1,2:9,10-Dibenzoheptazethrene (5), and 5,6:13,14-Dibenzoheptazethrene (6) Showing Quinoid Kekulé and Biradical Resonance Formsa

The benzene rings in red represent Clar’s aromatic sextet rings.

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 ). 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 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 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 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 represents the non-Kekulé biradical form of the Clar goblet.

Chart 3

Structures of p-Quinodimethane-Linked Bisphenalenyl 7 (Quinoid Kekulé and Biradical Resonance Forms) and Clar Goblet (8)a

The benzene ring in red represents Clar’s aromatic sextet ring.

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 –Chart ) 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 Kekulé or 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 –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 C2 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 two-electron 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 (NU)[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 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 11A state of structure 5 shows imaginary frequencies in C2 symmetry, which demonstrates that all other structures are minima. Structure 5 possesses two out-of-plane imaginary frequencies (b and a) in C2 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 C1 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 C1 symmetry would have been too costly. Since the difference between the planar C2 and nonplanar C1 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 C2 symmetry.

A wave function stability analysis[67] of the Kohn–Sham determinant[68] has been performed for the optimized 11A 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 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 (13B for structures 4-6 and 13B2 for structures 7 and 8) have also been optimized separately. In the spin-unrestricted Kohn–Sham formalism used, especially the low-spin BS solution 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 expressionwhere ET, ES, and EBS represent the energy of the triplet, singlet, and the BS solutions, respectively. ⟨ST2⟩ and ⟨SBS2⟩ 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 NU, we have also computed the multiple diradical character[70,71]y (i = 0, 1, 2) where 0 ≤ y ≤ 1 and y ≥ y from the spin-projected UHF (P-UHF) theory aswhere T is the orbital overlap between the corresponding orbital pairs which can be expressed by where nHONO– and nLUNO+ are the occupation numbers of the ith highest occupied NO (HONO) and the ith lowest unoccupied NO (LUNO) computed from UHF NOs. y0 = 0 indicates pure closed-shell and y0 = 1 indicates a pure diradical character. A perfect diradical has y0 = 1 and y1 = 0. Comparable y0 and y1 values indicate that in addition to the HONO/LUNO pair, non-HONO/LUNO pairs are also important. In the case of y0 = 1 and a large y1 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

Phenalenyl-Based Triangular Radicals

Figure shows the MO occupation schemes for the ground state of the phenalenyl-based neutral radicals. The symmetry is given both in C2 and D3 (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]

Figure 1

Molecular orbital (MO) occupation schemes for (a) the 12A2 (12A1″) state of phenalenyl (1), (b) the 13B2 (13A2′) state of triangulene (2), and (c) the 14B1 (14A2″) state of the π-extended triangulene system (3). The symmetry is shown as C2 (D3) notation.

In Table 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 C2 symmetry instead of the actual D3 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 2A1″ (D3 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.

Table 1

Excitation Energies ΔE (eV) between the Ground and the Degenerate First Excited States for the Phenalenyl-Based Neutral Radicalsa1–3 for Three Different Methods Using the 6-31G* and 6-311G(2d) Basis Setsb

system	orbital scheme	CASSCF	MR-CISD	MR-CISD+Q
Doublet–doublet Excitation (12A1″ → 12E″)
1	(σ)+π-spacec/ 6-31G*	3.190	3.076	2.864
	π-space/6-31G*	3.193	2.790	d
	π-space/6-311G(2d)	3.165	2.777	d
Triplet–Triplet Excitation (13A2′ → 13E′)
2	(σ)+π-spacec/ 6-31G*	2.702	2.696	2.709
	π-space/6-31G*	2.738	2.750	2.746
	π-space/6-311G(2d)	2.714	2.709	2.697
Triplet–Singlet Splitting (13A2′ → 11E′)
2	(σ)+π-spacec/ 6-31G*	0.842	0.798	0.703
	π-space/6-31G*	0.844	0.792	0.599
	π-space/6-311G(2d)	0.834	0.767	0.565
Quartet–Quartet Excitation (14A2″ → 14E″)
3	π-space/6-31G*	3.313	3.243	3.036
	π-space/6-311G(2d)	3.290	3.217	3.003
Quartet–Doublet Splitting (14A2″ → 12E″)
3	π-space/6-31G*	0.928	0.659	0.402
	π-space/6-311G(2d)	0.896	0.633	0.375

The symmetry is given as D3 notation, D3 (C2) 1: 2A1″ (2A2) → 2E″ (2B1,2A2); 2: 3A2′ (3B2) → 1,3E′ (1,3A1, 1,3B2); 3: 4A2″ (4B1) → 2,4E″ (2,4A2, 2,4B1).

For 1 and 3 a π-CAS(7,7)c and for 2 a π-CAS(4,7)c reference space are used, respectively.

Orbital occupation specifications are given in Table S1.

Contains intruder states.

The next higher non-Kekulé-type triangular PAH investigated is triangulene (2). Though it has an even number of carbon atoms, no Kekulé-type formula can be given. 2,6,10-Tri-tert-butyltriangulene[14] is the first example of an experimentally 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 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 . The 13E′ state is ∼2.7 eV higher than the 13A2′ 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 4E″ 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 ) 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 show that phenalenyl (1) has one (1a1″), triangulene (2) has two (degenerate 4e″), and π-extended triangulene (3) has three (one in 2a1″ and two in 6e″) singly occupied NO(s) (SONOs) in agreement with the MO occupation scheme given in Figure . 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).

Figure 2

Natural orbital occupation of the 12A1″ state of phenalenyl (1), the 13A2′ state of triangulene (2), and the 14A2″ state of the π-extended triangulene system (3) using a π-MR-AQCC/CAS(7,7)/6-311G(2d) for 1 and 3 and a π-MR-AQCC/CAS(4,7)/6-311G(2d) for 2.

In Figure 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 .

Figure 3

Density of unpaired electrons for the (a) 12A1″ state of phenalenyl (1) (NU = 1.33 e); (b) 13A2′ state of triangulene (2) (NU = 2.50 e), and (c) 14A2″ state of π-extended triangulene system (3) (NU = 3.66 e) using a π-MR-AQCC/CAS(7,7)/6-311G(2d) for 1 and 3 and a π-MR-AQCC/CAS(4,7)/6-311G(2d) calculation for (2) (isovalue 0.005 e bohr–3) with individual atomic population computed from Mulliken analysis.

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, NU presented in Figure , 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 non-negligible contributions generated from NOs other than the SONO(s) (Figure ).

Zethrenes, p-Quinodimethane-Linked Bisphenalenyl, and the Clar Goblet

Figures S1 and S2 show selected bond distances of the singlet ground state of structures 4–6 and 7–8 as computed at the UDFT/B3LYP/6-31G* level. The pronounced bond length alternation of the central p-quinodimethane subunit of the structures 4 and 5 indicates that the quinoid resonance form (Chart ) has a dominant 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 C2- and C1-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 (NU) 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 ) as shown in Figure S1 for structures 4 and 5 are very close to the experimentally observed[32] ones with a maximum deviation of 0.01 Å. The optimized structure of 7 (Figure S2) shows that the bond length labeled d (1.456 Å), which connects the two phenalenyl rings, and the bond lengths b (1.446 Å) and c (1.396 Å) of the central benzene ring are very close to the experimentally observed[35] bond lengths of 1.450, 1.437, and 1.394 Å, respectively.

Singlet–triplet splitting is investigated next in order to continue the discussion of polyradical character of structures 4–8. In Table , computed vertical and adiabatic (in parentheses) singlet–triplet splitting energies ΔE(S–T) are listed and compared with available experimental data.

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 ΔEproj 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) in Comparison with Available Experimental Dataa

system	ΔE(S – T)	MR-CISD	MR-CISD + Q	MR-AQCC	UDFT/B3LYP	⟨SBSi2⟩	ΔEproj	expt
4	E(13Bu) – E(11Ag)	1.188 (0.944)	0.904 (0.612)	0.927 (0.699)	0.598 (0.424)	0.228	0.672
5	E(13Bu) – E(11Ag)	0.776 (0.600)	0.516 (0.366)	0.671 (0.527)	0.384 (0.296)	0.602	0.540
6	E(13Bu) – E(11Ag)	0.354 (0.204)	0.447 (0.361)	0.414 (0.324)	0.201 (0.164)	0.889	0.351	0.165c
7	E(13B2u) – E(11Ag)	0.342 (0.244)	0.273 (0.153)	0.308 (0.197)	0.234 (0.211)	0.661	0.343	0.211d
8	E(13B2u) – E(11Ag)	0.026b (0.009)	0.075b (0.075)	0.170b (0.159)	0.038 (0.037)	1.159	0.085

In the MR calculations, π-CAS(4,4)b reference space and the 6-31G* basis set were used.

For 8, a RAS(1)/CAS(4,4)/AUX(1) reference space was used where RAS refers to a restricted active space and AUX to an auxiliary space (see the SI). Orbital occupation specifications are given in Table S3. In the multireference calculations, all σ orbitals were frozen.

Reference (32).

Reference (35).

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 ) 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 triplet state. These features are consistent with the theoretically observed small ΔE(S−T) of structure 6. Similar features are also observed for structure 7 where the adiabatic ΔE(S–T) computed at the same level is 0.197 eV, very close to the experimentally observed value[35] of 0.211 eV. For structure 8, the computed ΔE(S−T) is very small (mostly <0.1 eV, Table ) 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 together with the expectation values of S2 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 , ΔEproj). 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, ΔEproj of the nonplanar unrestricted C1 structure is 0.568 eV, very close that of the value of 0.540 eV calculated for the planar C2 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 . 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.

Figure 4

Natural orbital (NO) occupation of the 11Ag ground state of heptazethrene (4), 1,2:9,10-dibenzoheptazethrene (5), 5,6:13,14-dibenzoheptazethrene (6), p-quinodimethane-linked bisphenalenyl (7), and Clar goblet (8) obtained from π-MR-AQCC/CAS(4,4)/6-311G(2d) calculations for structures 4–7 and π-MR-AQCC/RAS(1)/CAS(4,4)/AUX(1)/6-311G(2d) calculation for structure 8.

Comparison of the HONO/LUNO occupation numbers (Table ) 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 closed-shell 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 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.

Table 3

Comparison of the NO Occupation for the Singlet State of the Structures 4–8 Obtained from π-MR-AQCC and UHF Calculations, Respectively

	HONO		LUNO
system	π-MR-AQCC	UHF	π-MR-AQCC	UHF
4a	1.705	1.234	0.295	0.766
5a	1.624	1.173	0.375	0.827
6a	1.387	1.118	0.610	0.882
7a	1.450	1.116	0.558	0.883
8b	1.080	1.015	0.917	0.985

π-MR-AQCC/CAS(4,4)/6-311G(2d) and UHF/631G* calculations.

π-MR-AQCC/RAS(1)/CAS(4,4)/AUX(1)/6-311G(2d) and UHF/6-31G* calculations.

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 ). 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 and Figure 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 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 . 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 and 3.

Figure 5

Density of unpaired electrons for the 11Ag state of (a) heptazethrene (4) (NU = 1.03 e), (b) 1,2:9,10-dibenzoheptazethrene (5) (NU = 1.49 e), and (c) 5,6:13,14-dibenzoheptazethrene (6) (NU = 2.24 e) using the π-MR-AQCC/CAS(4,4)/6-311G(2d) approach (isovalue 0.003 e bohr–3) with individual atomic population computed from Mulliken analysis.

Figure 6

Density of unpaired electrons for the 11A state of (a) p-quinodimethane-linked bisphenalenyl (7) (NU = 1.86 e) using the π-MR-AQCC/CAS(4,4)/6-311G(2d) method and (b) Clar goblet (8) (NU = 2.88 e) using the π-MR-AQCC/RAS(1)/CAS(4,4)/AUX(1)/6-311G(2d) method (isovalue 0.003 e bohr–3) with individual atomic population computed from Mulliken analysis.

The unpaired density for structure 7 (Figure ) shows a pattern which is more delocalized than the one indicated by the two mesomeric forms (7a′ and 7b′) given in Chart . Thus, in this case, the unpaired density can be better represented 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 ), the unpaired character is mostly located at the zigzag edges with the largest contribution at their centers (position at C9/18, see Chart ). 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 compares the NU values computed at π-MR-AQCC level with the multiple diradical character indices, y (i = 0,1,2..), obtained from P-UHF theory for structures 4–8. To obtain a common basis for comparison with the NU’s, the y values were multiplied by a factor of 2. It is observed that for the singlet ground states of structures 4 and 5 the 2·y0 values, which are computed from the HONO/LUNO UHF occupations, are almost twice of the NU 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 ). The y0 values reported in ref (32) for the structures 4–6 are somewhat smaller as compared to our values. However,, the trend of increasing y0 value as one moves from structure 4 to 6 is the same. For structure 5, the y0 value as computed for C1 optimized structure is 0.649, very close to that of the C2 planar structure of 0.664. This indicates the similarity of the NO occupation numbers between the two structures. Comparing the NU values derived from different NO selections, it is noted that the total NU 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 ), 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 (TCNE22–) and neutral K2TCNE2 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 NU 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 NU 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

Comparison of the Number of Effectively Unpaired Electrons (NU) with the Multiple Diradical Characters 2y0 and 2y1 for Singlet and Triplet States of Structures 4–8

	singlet state						triplet statea
	π-MR-AQCCb				P-UHFc		π-MR-AQCCb
system	NU	NU from HONO/LUNO	NU from HONO-1/LUNO+1	NU from non- HONO/LUNO	2y0	2y1	NU	NU from HONO/LUNO	NU from non- HONO/LUNO
4	1.026	0.506	0.124	0.520	1.112	0.182	2.589	1.990	0.599
5	1.492	0.744	0.185	0.748	1.328	0.357	2.795	1.996	0.799
6	2.241	1.442	0.200	0.799	1.534	0.262	2.750	1.998	0.752
7	1.865	1.283	0.107	0.582	1.540	0.050	2.598	1.999	0.599
8	2.880	1.974	0.257	0.906	1.940	0.410	2.816	2.000	0.816

2y0 for the triplet state is two.

The 6-311G(2d) basis and a π-CAS(4,4) reference space for 4–7 and RAS(1)/π-CAS(4,4)/AUX(1) for 8 were used.

P-UHF with 6-31G* basis set was used.

Table shows that for all structures the triplet state maintains practically a constant NU value (from 2.589 e to 2.816 e), whereas for the singlet state of structures 4–8 a large change in the NU 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 NU 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 and 7). These differences are naturally larger for the cases with smaller NU 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 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 NU 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 and Figure S3) the distribution of unpaired density is very similar in nature.

Figure 7

Density of unpaired electrons for the 13Bu state of (a) heptazethrene (4) (NU = 2.59 e); (b) 1,2:9,10-dibenzoheptazethrene (5) (NU = 2.80 e), and (c) 5,6:13,14-dibenzoheptazethrene (6) (NU = 2.75 e) using the π-MR-AQCC/CAS(4,4)/6-311G(2d) approach (isovalue 0.003 e bohr–3) with individual atomic population computed from Mulliken analysis.

Conclusions

In this work two different types of PAHs have been studied (a) non-Kekulé triangular structures with a high-spin ground state and (b) PAHs with singlet polyradical character. For the first case, phenalenyl (1), triangulene (2), and a π-extended triangulene system (3) have been chosen. In the second case, a series of three zethrenes, heptazethrene (4), 1,2:9,10-dibenzoheptazethrene (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 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. Spin-projection 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 y0 and y1 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.