Can a Topological Approach Predict Spin-Symmetry Breaking in

Article pubs.acs.org/JPCA

Can a Topological Approach Predict Spin-Symmetry Breaking in Conjugated Hydrocarbons? Jean-Paul Malrieu* and Georges Trinquier Laboratoire de Chimie et Physique Quantiques, IRSAMC-CNRS-UMR5626, Université Paul-Sabatier (Toulouse III), 31062 Toulouse Cedex 4, France S Supporting Information *

ABSTRACT: The closed-shell mean-field single determinants of large alternant hydrocarbons are frequently unstable with respect to a possible spin-symmetry breaking which produces different orbitals for the α and β electrons, either in Hartree−Fock or in Kohn− Sham DFT calculations. The present work shows that one may easily predict whether such a symmetry breaking will take place from the elementary topological Hückel Hamiltonian which introduces a simple hopping integral t. The demonstration makes use of the simplest representation of the bielectronic repulsion, namely, the Hubbard bielectronic operator, reduced to an on-site repulsion U, and takes benefit of the mirror theorem. A recipe is proposed to determine the relevant t/U ratio for a given exchange-correlation potential. The symmetry-breaking phenomenon first concerns the mixing between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), but it may eventually run on other pairs of mirror orbitals. These symmetry breakings may take place while the other molecular orbitals keep a closed-shell character. The spin polarization of these MOs, appearing in typical unrestricted mean-field calculations, is an induced and amplifying effect, which has to be distinguished from the symmetry breaking itself. Special attention is paid to the possible appearance of multiple symmetry breakings, leading to a polyradical character. The model is tested on six series of polycyclic hydrocarbons. This elementary approach sheds new arguments on the debate concerning the dior polyradical character of polyacenes. place when the second factor prevails over the first one. Notice that the spin-symmetry breaking takes place in densityfunctional-theory calculations (DFT) as well as in Hartree− Fock descriptions (HF), although at different interatomic distances.17 Since the exchange-correlation potential reduces to some extent the electron−electron repulsion, the symmetry breaking takes place at larger interatomic distance in DFT than in HF, as will be illustrated in section 3.1. This is consistent with the fact that a treatment of correlation may cancel the symmetry breaking in the orbital optimization,18 In both DFT and HF treatments the working physical effect inducing this phenomenon starts with a HOMO/LUMO mixing. In complex systems, for instance in multiple bonds,19 but also in extended systems,20−22 the symmetry breakings may be multiple and of various characters (some of them may keep the desired spin properties, being S2 eigenfunctions, but breaking the desired space symmetry). The symmetry breaking also affects extended systems, in particular half-filled band systems.23,24 The present work concentrates on a specific type of such systems, namely, the conjugated hydrocarbons, where the π electrons in the valence π orbitals define a halffilled band problem. We only consider systems with 2n sites and alternant hydrocarbons. The alternant hydrocarbons accept a bipartition of the sites, attributing one of two colors, say red

1. INTRODUCTION The discussion of the occurrence of space and/or spin symmetry breakings of mean-field descriptions of electronic systems has a long story.1−8 Their appearance is clearly connected to the electron−electron correlation problem. The symmetry-adapted mean-field wave function may impose some equality between components of the wave function which intrinsically have very different energies. This is the case for the paradigmatic problem of the H2 molecule, which manifests a symmetry breaking of the Hartree−Fock single determinant at large enough interatomic distances. The symmetry imposes the coefficients of neutral and ionic valence-bond components of the closed-shell single determinant to be equal, while the ionic one is of much higher energy. The weight of this component should tend to zero when the interatomic distance increases. Breaking the symmetry of two spin orbitals enables one to diminish the weight of the ionic component but tends to localize the α-spin electron on one of the two atoms and the βspin electron on the other atom.9,10 This solution tends to take into account an important correlation correction, namely, the so-called left−right nondynamical correlation effect in this particular case. Of course the resulting wave function is no longer a pure singlet state, but a mixture of singlet and triplet states. This is the famous spin-contamination problem of broken-symmetry solutions, which may receive various more or less sophisticated corrections.11−16 The occurrence of symmetry breaking is governed by a competition between interatomic delocalization and electron−electron repulsion, and it takes © XXXX American Chemical Society

Received: July 28, 2016 Revised: November 3, 2016 Published: November 8, 2016 A

DOI: 10.1021/acs.jpca.6b07597 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A and blue, to each of the sites in such manner that the first neighbors of each atom of one color are of the other color. This simply implies that the molecular graph is free from oddmembered rings. Ovchinnikov has formulated a simple statement concerning the preferred spin multiplicity of the ground state, saying that this spin multiplicity is equal the difference between the numbers of red and blue sites, augmented by 1.25 If the two numbers are equal, the ground state is therefore expected to be a singlet state, and the present study will be restricted to this case actually. Ovchinnikov’s statement was established from a magnetic (Heisenberg) picture of the π electron population, that is, from the strong correlation limit, but it may be reached whatever the electronic correlation versus delocalization ratio.26,27 The magnetic properties of conjugated hydrocarbons attract an intense interest,28 especially in view of spintronics applications. The ferromagnetic architectures receive a special interest,28,29 but singlet diradicals also deserve attention from theoreticians30−36 and experimentalists.37,38 Spin-symmetry breakings are observed in large enough conjugated hydrocarbons, which have the same number of atoms of both colors and should be considered as of singlet ground state according to Ovchinnikov’s rule. The occurrence of symmetry breaking in these systems suggests an antiferromagnetic character, with a low-lying triplet state. Among them, some are clearly disjoint diradicals when it is impossible to proceed to a Lewis-type pairing of electrons in disjoint double bonds, that is, to write a Kekulé formula.39−41 Such systems require intrinsically openshell descriptions, and they will be discarded from the present work as we concentrate only on “Kekuléisable” graphs, usually expected to be of closed-shell character. In HF calculations, symmetry breaking appears in many polycyclic hydrocarbons as small as anthracene,17 (and even in benzene and naphthalene as will be discussed in section 3.3), while in the Kohn−Sham DFT paradigm (KS-DFT), it appears only for larger ones, as it will be seen below in the same section. At least two series of such systems are known to exhibit a symmetry breaking. The first one is the series of para-dimethylene-para-polyphenylenes, pR2C−(C6H4)n−CR2 and related architectures.42 This family is better known after the name of the chemists who first synthetized derivatives of this type, namely, Thiele (n = 1),43 Chichibabin (n = 2),44 and Müller (n = 3).45,46 The appearance of unpaired electrons on terminal CH2 groups is not surprising here since it maintains the aromaticity of the benzene rings, which would be destroyed by the full electron pairing on double bonds into a paraquinonic configuration. The appearance of instabilities in the series of linear polyacenes was more surprising and has been the matter of ongoing debates,47−51 as briefly recalled in section 3.2. The present work first establishes a distinction between the occurrence of a spin-symmetry breaking and the subsequent spin-polarization phenomenon. The symmetry breaking may concern one MO only, the other MOs keeping a closed-shell character. Once the symmetry breaking of a given MO has taken place, one may lower the mean-field energy by giving different space parts to the previously closed-shell MOs as well, introducing a spin-polarization phenomenon, which is a secondary effect, induced by the symmetry breaking. The present work makes use of very primitive tools, namely, the topological purely monoelectronic Hückel Hamiltonian,52 and the simplest representation of the electron−electron repulsion, provided by the Hubbard Hamiltonian,53 as being reduced to on-site interactions. The subsequent rationalizations involve

only two parameters, namely, the hopping integral t between nearest neighbors and the on-site electron−electron repulsion U. We shall show that using the Hückel MOs, taking benefit of the mirror theorem, one may predict whether a conjugated hydrocarbon may present a symmetry breaking, resulting from a mixing between the HOMO and the LUMO. The symmetrybreaking condition is expressed as a topology-dependent prefactor multiplied by a t/U ratio which only depends on the functional. A method to fix the value of this ratio from elementary computations on the ethylene molecule will be proposed in section 3.1. Systematic comparisons between our predictions and the behavior of UDFT-B3LYP calculations will be presented for six series of compounds in section 3.2, confirming the validity of our model. Moreover the method enables one to predict whether the symmetry breaking may also concern deeper occupied MOs, for instance the HOMO-1 and the LUMO+1 orbitals, as has been observed in very long polyacenes and in regular fused rhomboidal polyacenes.54 Some comparative UHF results will be presented in section 3.3. Expressed in a configuration-interaction approach (CI), the analysis makes it possible as well to estimate the correlation energy brought by double excitations from one occupied orbital to its mirror antibonding counterpart, whether this excitation generates a symmetry breaking or not. It furnishes an estimate of the occupation numbers of occupied and virtual orbitals. Of course, the appearance of symmetry breaking in mean-field calculations is frequently accepted as a proof for diradical character of the system. But this symmetry breaking may be of variable strength, the mean value of the spin operator S2 in the broken-symmetry single determinant may vary between zero (for a pure singlet state) to one (equal mixture of triplet and singlet sates for a diradical), and to higher values when other shells tend to get open. A quantitative measure of diradical character can be found in ref 13. Occupation numbers of occupied and virtual orbitals calculated from limited CI treatment offer a more reliable evaluation of the open-shell character of the singlet states. Simple topological transcription of spin-polarization effects in terms of configuration interaction avoids the difficulties of spin decontamination techniques based on ⟨S2⟩ values.11−15

2. FORMAL DEVELOPMENT 2.1. Simplified Tools for the Treatment of Conjugated Hydrocarbons. The simplest topological Hückel Hamiltonian only involves one hopping integral t (of negative sign) between adjacent sites.52 It assumes the energies of the (supposedly orthogonal) 2pz atomic orbitals to be the same for the conjugated carbons, and the CC bonds to have similar lengths (which is reasonable for fused polycyclic hydrocarbons). It takes the form Huc =

∑ t(ap+aq+aq+ap) p,q

The Hubbard Hamiltonian simply adds the on-site bielectronic repulsion through a positive parameter U Hub =

∑ t(ap+aq+aq+ap) + U ∑ np↑np↓ p,q

p

which only acts when two electrons occupy the same atomic orbital.53 We shall not search the self-consistent single determinant solution of the Hubbard Hamiltonian, since it requires, in principle, a black-box computation. We shall use the B


Article

The Journal of Physical Chemistry A Hückel orbitals to build a closed-shell reference determinant, occupying the n lowest-energy MOs in a 2n-site problem. The Hückel orbitals are linear combinations of the atomic orbitals p φi =

that is, they are solutions of the Hartree−Fock equations for the Hubbard Hamiltonian. 2.2. The Symmetry-Breaking Condition. Φ0 only interacts with doubly excited determinants. One may concentrate first on the doubly excited determinant obtained by exciting the two electrons of the highest-occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO). But for sake of broader generality we consider the double excitations from an MO φi to its mirror MO φi*. The interaction between Φ0 and the closed-shell doubly excited determinant is

∑ cipp p

Huc|φi⟩ = εi|φi⟩

and have an energy εi. The closed-shell description of the ground state is

Φ0 = |Πiφφ | i i̅

⟨Φ0|H |al+*a l+̅ *a i ̅ ai Φ0⟩ = K i * i

It is now important to recall the “Mirror Theorem”, established by Longuett-Higgins55 and Hall,56 which evidence a two-by-two correspondence between occupied and virtual unoccupied orbitals of the Hückel Hamiltonian. The first part of the theorem establishes a correspondence between each occupied MO |φi⟩ of energy εi and a virtual MO |φi*⟩, the energy of which is −εi.

which can be expressed as ⟨Φ0|H |al+*a l+̅ *a i ̅ ai Φ0⟩ = U

Jii = U

cip4

This interaction is positive and expected to be larger than for the more numerous four-open shell excitations ⟨Φ0|H |al+*ak+*a j ̅ ai Φ0⟩ = U

q′

cipcjpcl * pck * p

where the products of coefficients may have random signs. Actually, the largest interaction should concern the double excitation from the occupied MO presenting largest coefficients, that is, the largest number of nodes, hence the singly occupied molecular orbital (SOMO), to its mirror orbital, the LUMO. The energy difference between Φ0 and the closed-shell doubly excited determinant al+*a l+̅ *a i ̅ ai Φ0 is given by

⟨Φ0|H |aj+*ai Φ0⟩ = ⟨φi|F |φj *⟩

⟨al+*a l+̅ *a i ̅ ai Φ0|H |al+*a l+̅ *a i ̅ ai Φ0⟩ − ⟨Φ0|H |Φ0⟩

with F = h + ∑k = 1, n 2Jk − Kk , h being the monoelectronic part of the Hamiltonian (i.e., the Hückel Hamiltonian), J and K being the Coulomb and exchange operators, respectively. Hence ⟨Φ0|H |aj+*ai Φ0⟩ = ⟨φi|Huc|φj *⟩ + U

∑

= 2(εi * − εi) + Jii + Ji * i * − 4Jll * + 2K ii *

It is easy to see that, due to the mirror theorem, the bielectronic terms cancel for the Hubbard Hamiltonian, so that

cipcj * p

⟨al+*a l+̅ *a i ̅ ai Φ0|H |al+*a l+̅ *a i ̅ ai Φ0⟩ − ⟨Φ0|H |Φ0⟩ = 2(εi * − εi)

p = 1,2n

∑

∑ p = 1,2n

when going from an occupied MO to its mirror empty MO, the sign of the coefficients changed for the atoms of one of the two colors, remaining the same for the other-color sites. Notice that the determinant Φ0 does not interact with the singly excited determinants, obtained by sending one electron from an occupied orbital to an antibonding one

= −4εi

ckp 2)

k = 1, n

If one concentrates on the interaction between these two determinants, one defines a 2 × 2 CI matrix

The monoelectronic part is zero since the orbitals are eigenfunctions of the Hückel operator. For a half-filling situation (2n electrons in 2n orbitals) the atomic charges qp = 2(∑k = 1, n ckp 2) of the Hückel solution are equal to 1, again

Φ0

E0

Jii

al+*a l+̅ *a i ̅ ai Φ0 Jii E0 − 4εi

due to the mirror theorem. Since the quantity

∑

∑

⟨Φ0|H |al+*a l+̅ *a i ̅ ai Φ0⟩ = Jii

∑ cipp − ∑ ciq′q′

× (2

p = 1,2n

hence

q′

p

cip4

p = 1,2n

∑ cipp + ∑ ciq′q′

φi* =

∑

As Coulomb repulsion of the two electrons in the MO φi has the same value,

The second part of the theorem concerns the coefficients of these orbitals. If one divides the sites (and therefore the AOs) into the red and the blue subsets, giving a prime symbol to the second set only, the bonding and antibonding mirror MOs may be written p

cip2ci2* p = U

p = 1,2n

Huc|φi*⟩ = −εi|φi*⟩

φi =

∑

This matrix is reminiscent of the elementary problem of H2 in the minimal basis set, the MOs φi and φi* playing the role of the MOs σg and σu, respectively. The symmetry-breaking problem in the two-electron bond is a textbook question. One may actually introduce the “localized” MOs

cipcj * p = ⟨φi|φj *⟩ = 0

p = 1,2n

the Hückel MOs satisfy the Brillouin’s theorem ⟨Φ0|H |aj+*ai Φ0⟩ = 0

φai = (φi + φi *)/ 2 C


Article

The Journal of Physical Chemistry A

mirror antibonding MO may induce a symmetry breaking in the mean-field calculation of ground-state determinant. Of course, this criterion has to be applied first to the HOMO/ LUMO pair since the corresponding orbital energy is the lowest in absolute value. But, as will be shown below, in some highly delocalized hydrocarbons, a second symmetry breaking may take place involving the HOMO-1 and LUMO+1 pairs, or even deeper pairs. It may be interesting as well to recall14 that symmetry breaking stabilizes energy by the quantity

φbi = (φi − φi *)/ 2

which, according to the mirror theorem, are respectively defined on red and blue color sites. The limited-CI problem may again be formulated in terms of interaction between neutral configuration ΦN (one electron in each of these two “localized” orbitals), |ΦNi⟩ = |Πj ≠ iφφ (φ φ̅ + φbiφai̅ )/ 2 | j j̅ ai bi

and ionic configuration ΦI (two electrons in one of these two “localized” orbitals),

ΔE BS = −

|ΦIi⟩ = |Πj ≠ iφφ (φ φ + φbiφbi̅ )/ 2 | j j̅ ai ai̅

The interaction between these two valence-bond type configurations is

and that the mean value of the S operator relative to the broken-symmetry solution is ⟨S2⟩ = 1 −

and their energy difference is Ui = Jai , ai − Jai , bi

εi 2 Jii 2

(4)

4εi 2 + Jii 2

(5)

In this way, one may have an upper bound to correlation energy by summing all independent-pair excitation corrections, and one may have an upper bound of the occupation number of the corresponding occupied orbital as

It is known that one may break the symmetry of a given MO introducing different space parts for the α and β spin electrons, here (1)

2

ni =

φi̅″ = sin θφi̅ + cos θφi̅ * = sin αφai + cos αφbi

1+

in a “frozen core” broken-symmetry single determinant,

2

( ) δii → i * i * Jii

Of course, double excitations between nonmirror orbitals ii → j*j*, not inducing intrinsic symmetry-breaking, contribute to both correlation energy and lowering the occupation numbers of occupied orbital. The instability conditions concerning the possible mixing by an angle φ of the HOMO and LUMO might be extended to nonalternant hydrocarbons. They imply a number of quantities, namely, the energies of the HOMO h (εh) and of the LUMO l (εl), the Coulomb integrals Jhh, Jll, and Jhl, the exchange integral Klh, and the integrals ⟨h|Jl|l⟩ and ⟨h|Jl|l⟩. All of these quantities are easily expressed from the Hubbard Hamiltonian. Taking the energy of the closed-shell determinant is zero of energy, and introducing the excitation energies

′φ″| |ΦBS, n⟩ = |Πj ≠ iφφφ j j̅ i i̅ |ΦBS, n⟩ = |Πj ≠ iφφ ((cos α)2 φaiφbi̅ + cos α sin α j j̅ × (φaiφai̅ + φbiφbi̅ ) + (sin α)2 φbiφai̅ )|

One sees immediately the qualitative change introduced by the symmetry-breaking; it reduces the weight of the ionic VB components of the wave function, to the price of a mixing with a triplet component since the coefficients of the two neutral components are no longer equal. In this determinant, the symmetry breaking is restricted to a single electron pair. One may call this phenomenon an intrinsic symmetry breaking, taking place before the spin polarization of the other MOs. This distinction has already been underlined in the literature.57−59 It has been recently identified and used by Coulaud et al. to discriminate between the so-called kinetic exchange and the spin-polarization mechanism.57,58 The condition for the appearance of the symmetry breaking is well-known (see for instance ref 14); it takes place provided that the ratio

Δhh → ll = 2(εl − εh) + Jll − Jhh Δh → l = εl − εh + Jhl − Jhh + 2Khl

the energy of the broken-symmetry determinant becomes E(φ) = sin φ4Δhh → ll + (sin 2φ 2Δh → l + sin 2φ(cos φ 2⟨h|Jh |l⟩ + sin φ 2⟨h|Jl |l⟩)/2

< 1/2

which can be derived with respect to the angle φ. 2.3. Spin Polarization of Inner Shells. The above derivations treat the symmetry breaking as concerning essentially two orbitals, as actually implicitly assumed in the famous Yamaguchi’s correction.11−13 The mean value of S2 could not exceed one if two MOs only were concerned. In practice, the symmetry breaking involving mainly the HOMO and the LUMO induces a spin polarization of the other orbitals,

that is, in the Hubbard model, provided that the ratio

|εi| 2 occupied orbitals. It would proceed through a localization of the m highest canonical orbitals into a set of m localized orbitals φ″k . The transformation should maximize the repulsion between the electrons,

∑

Figure 1. Representative examples in the families of conjugated hydrocarbons presently addressed. From top to bottom: p-dimethylene p-polyphenylene (size [3]), periacene (size [4], peritetracene), benzo-periacene (size [5], benzo-peritetracene), rhombus (size [4]), ”spiked” rhombus (size [4]), polyacene (size [6], hexacene).

DFT version of the Gaussian 09 quantum chemistry package, with standard B3LYP hybrid functional and 6-311G** basis set.61,62 Some unrestricted Hartree−Fock (UHF) calculations have been reported in Section 3.3 for comparison, using the same basis set and the same code. The ratio |U/t| plays a key role and must be fixed. It has no reason to be equal in HF and DFT calculations since the exchange-correlation potential is supposed to introduce some correlation effects which are absent in HF. Moreover the appropriate |U/t| ratio is functional-dependent. We suggest here a simple procedure to fix its value. It consists in two independent calculations on the ethylene molecule. The difference between the energy of the bonding π MO, Fππ, and that of the antibonding π* MO, Fπ*π* may be exploited to fix the t value

−1 ⟨φk″φk″|r12 |φk″φk″⟩max

k = 1, m

One recognizes here the well-known criterion proposed by Ruedenberg to localize the SCF occupied MOs.60 Instead of applying it on all of the occupied orbitals, which would lead to bond orbitals, the energies of which are very large, and which cannot generate symmetry breaking, we recommend to apply the localization to a set of orbitals close to the Fermi level. Actually, the appearance of several instabilities is ruled by a conflict between orbital energies, which must remain small, hence involving only high-lying orbitals, and the repulsion of the electron pair occupying the same orbital. Again, these rotations are relevant only when they lead to localized orbitals.

2t = Fππ − Fπ * π * In order to treat the two MOs on the same level, these energies are extracted from the calculation of the 3Bu π→π* excited state, described as a single determinant 3Φ = |Φσππ*|, where Φσ represents the mean field description of the σ core. Going now to the closed-shell description Φ0 = |Φσππ̅| of the 1Ag ground state, the energy difference between its energy and that of the 3 Bu state is easily calculated in the Hubbard model as

3. NUMERICAL TESTS 3.1. Computational Details and Determination of a Relevant |U/t| Ratio. The predictions of our topological model are now put to test in typical families of polycyclic conjugated hydrocarbons (Figure 1), which have been the object of UDFT calculations. The results are obtained from the G


Article

The Journal of Physical Chemistry A Table 2. Symmetry Breaking in Periacenesa

E(Φ0) − E(2Φ) = 2t + U /2

These two equations fix the values of the two parameters t and U for a given geometry of the ethylene molecule, and one may propose a geometry-dependent Hubbard Hamiltonian.63 Since the model Hamiltonians calculation fix an ideal geometry and since most of the hereafter considered systems are polycyclic hydrocarbons, with a mean value of the C−C bond length close to 1.40 Å, this distance will be accepted. The B3LYP calculation leads to t = −0.103 au and U = 0.125 au, that is, a ratio U/|t| = 1.23, while the corresponding quantities are t = −0.126 au and U = 0.294 au, that is, a ratio U/|t| = 2.33 in the HF treatment. This ratio is much larger in HF than in DFT, as expected, and as systematically verified in the study of antiferromagnetic dimers.64 If one twists the CH2 groups by a dihedral angle θ, the RDFT solution is stable up to θ = 52°. At this angle the hopping integral is reduced to the value t′ = t cos θ, the ratio U/t′ = U/(t cos θ) = 2.02, which is almost the critical value (2.0) for the appearance of the symmetry breaking. If one accepts that the DFT calculation incorporates the dynamical correlation effects, it is possible to compare its predictions with those extracted from sophisticated ab initio configurationinteraction calculations, treating the impact of the dynamical correlation on energy differences.63,65 For the 1.40-Å CC bond distance, ref 63 reports the values t = −0.101 au and U = 0.127 au, in surprising agreement with the B3LYP DFT estimates. The resulting U/|t| ratio (1.26) is extremely close to the abovereported UDFT estimate. Notice that in the above case the value of U is the difference between the on-site repulsion Jaa and the repulsion between adjacent atoms Jab and U = Jaa − Jab. A somewhat larger value of U should be used in our modeling, since Jab, which concerns sites of opposite colors, does not play a role in the evaluation of Jii. We have therefore adopted a somewhat larger value U/|t| = 1.5 for the comparisons between the predictions of our model and the results of the B3LYP computations in our table. 3.2. Numerical Tests with the B3LYP ExchangeCorrelation Functional. Our predictive model has been tested on six series of polycyclic conjugated hydrocarbons. Tables 1−6 report on the left columns the characteristic values

n

orbitals

−εi

Jii

−εi/Jii

⟨S2⟩

4

HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR

0.095 0.493 0.294 0.053 0.347 0.200 0.032 0.247 0.139 0.020 0.117 0.068

0.114 0.070 0.124 0.113 0.062 0.118 0.109 0.059 0.115 0.104 0.058 0.113

0.83 7.09 2.37 0.47 5.55 1.69 0.29 4.16 1.21 0.19 3.07 0.87

0.99

5

6

7

1.16

1.28

1.53

Energies are in |t| units (assuming U = 1.5|t|); LR refers to left−right localized combinations of HOMO and HOMO-1; bold numbers highlight critical ratios less than 1, suggesting shell opening; the last column refers to UDFT calculations.

a

Table 3. Symmetry Breaking in Benzoperiacenesa n

orbitals

−εi

Jii

−εi/Jii

⟨S2⟩

4 5 6

HOMO HOMO HOMO HOMO-1 LR HOMO HOMO-1 LR

0.150 0.091 0.057 0.292 0.174 0.037 0.216 0.126

0.108 0.109 0.108 0.061 0.110 0.106 0.058 0.109

1.39 0.83 0.53 5.00 1.58 0.35 3.80 1.16

0. 0.90 1.14

7

a

1.30

Same comments as in Table 2.

Table 4. Symmetry Breaking in Rhombuses of Sizes [n] modela n

orbitals

−εi

Jii

−εi/Jii

⟨S ⟩

Δα

3 4

HOMO HOMO HOMO-1 HOMO HOMO-1 LR

0.186 0.071 0.271 0.025 0.125 0.075

0.118 0.106 0.080 0.102 0.070 0.054

1.58 0.67 3.39 0.25 1.79 0.72

0. 0.79

0.051 0.027

1.50

0.013

5

Table 1. Symmetry Breaking in p-Dimethylene Polyphenylenes H2C−(C6H4)n−CH2a

UDFTb 2

Same comments as in Table 2. b⟨S2⟩ refers to the ms = 0 solution; Δε corresponds to (α2 − α1), the energy gap (in au) between highest α occupied levels for the ms = 1 solution.

a

n

−εi

Jii

−εi/Jii

⟨S ⟩

1 2 3 4 5

0.311 0.135 0.064 0.031 0.016

0.361 0.272 0.242 0.230 0.226

0.86 0.50 0.26 0.13 0.07

0.b 0.32 0.98 1.07 1.08

2

Table 5. Symmetry Breaking in Spiked Rhombuses of Sizes [n]a

Energies are in |t| units (assuming U = 1.5|t|); last column refers to UDFT calculations. bSymmetry breaking observed only at triplet geometry.

a

n

orbitals

−εi

Jii

−εi/Jii

⟨S2⟩

2

HOMO HOMO-1 HOMO HOMO-1 LR HOMO HOMO-1

0.097 0.445 0.029 0.271 0.107 0.008 0.071

0.296 0.160 0.283 0.080 0.174 0.280 0.106

0.33 2.77 0.10 3.39 0.61 0.03 0.71

0.80

3

of our topological model, that is, the energies εi of the highest occupied MOs (in t units), the Coulomb integrals Jii of these MOs in −t units, assuming U = −1.5t, the ratio of these quantities, and in the right columns the numerical results of the UDFT calculations relative to the ms = 0 solution, namely, the mean value of S2 for the lowest-energy solution and in Tables 4 and 6 the HOMO−LUMO gap in au. These calculations concern optimized geometries relative to the ms = 0 solution, which are available upon request.

4 a

1.30

2.06

Same comments as in Table 2.

p-Dimethylene Polyphenylenes. In this series, also named p-dimethylene quinodimethanes, p-CH2−(C6H4)n−CH2, the external CH2 groups may either enter in CC double bond, to the price of a complete quinonization of the intermediate H


Article

The Journal of Physical Chemistry A Table 6. Symmetry Breaking in Polyacenes of Sizes [n] modela

UDFTb

n

orbitals

−εi

Jii

−εi/Jii

⟨S ⟩

Δε (ms = 1)

Δε (ms = 2)

7 8

HOMO HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR HOMO HOMO-1 LR

0.134 0.109 0.347 0.228 0.090 0.295 0.192 0.075 0.253 0.164 0.064 0.220 0.142 0.055 0.192 0.124 0.048 0.169 0.108 0.042 0.150 0.096 0.037 0.134 0.086

0.114 0.104 0.084 0.140 0.097 0.078 0.130 0.090 0.074 0.121 0.084 0.070 0.114 0.078 0.066 0.108 0.073 0.063 0.101 0.070 0.059 0.096 0.066 0.057 0.091

1.17 1.05 4.13 1.63 0.93 3.78 1.48 0.83 3.42 1.36 0.76 3.14 1.25 0.71 2.91 1.15 0.66 2.68 1.07 0.60 2.54 1.00 0.56 2.35 0.94

0. 1.05

0.021 0.014

0.090 0.076

1.24

0.009

0.064

1.39

0.005

0.055

1.54

0.002

0.047

1.71

0.001

0.040

1.90

0.003

0.034

2.10

0.005

0.029

2.30

0.006

0.024

9

10

11

12

13

14

15

2

Same comments as in Table 2. b⟨S2⟩ refers to the ms = 0 solution; Δε corresponds to energy gaps (in au) between highest α occupied levels, namely, (α2 − α1) for ms = 1 solution, and (α4 − α1) for the ms = 2 solution. a

benzene rings, or keep an unpaired electron and leave the benefit of aromaticity to these rings. The first member (n = 1) present a closed-shell singlet in DFT calculations, at least after geometry optimization. Our model predicts a symmetry breaking (Table 1). The contradiction is not so severe, since at the geometry of the triplet state the ms = 0 is symmetrybroken.42 Our model makes use of equal values of the hopping integrals, that is, equal bond lengths, and this is appropriate for the triplet state geometry, not for the quinonic geometry which presents strong contrasts between single and double bonds. The next members are subject to symmetry breaking with increasing ⟨S2⟩ values, in agreement with the model, which diagnoses a single symmetry breaking by HOMO/LUMO mixing. The small excess of ⟨S2⟩ beyond 1 is due to the spin polarization. Periacene Nanoribbons. Although it may have escaped to previous works,66 peritetracene (or tetra-benzo-coronene) exhibits a symmetry breaking in the above-specified DFT calculations, with ⟨S2⟩ = 0.99. This agrees with our model, which predicts instability with respect to HOMO/LUMO mixing, since the critical ratio, |εi|/Jii = 0.83, is smaller than 1 (Table 2). Referring to Clar’s theory of aromatic systems,67 and some of its recent applications,68−75 one understands this molecule may exhibit five sextets, to the price of creating two allyl-type radicals at the top and bottom centers (Scheme 1). More generally, this molecule belongs to the periacene series,76−86 that is, nanoribbons of width three, with n rings in the upper and lower ranks and n − 1 rings in the central rank. As shown in Table 2, a second symmetry breaking may take place for n = 7, but it is not a direct instability with respect to the HOMO-1/LUMO+1 mixing. It is obtained only after a

Scheme 1. Maximizing the Number of Clar’s Sextets in Peritetracene

left/right localization of these orbitals, according to the procedure defined above in Section 3.3. This second symmetry breaking would suggest to view periheptacene as in Scheme 2, and actually UDFT calculations suggest this system is an intermediate between a diradical and a tetraradical, since the Scheme 2. Maximizing the Number of Clar’s Sextets in Periheptacene

I


Article

The Journal of Physical Chemistry A value of ⟨S2⟩ is calculated at 1.53. However, the ms = 2 solution is here lying at 14 kcal/mol above the ms = 0 solution. In view of this large energy difference, the di- vs tetraradical character remains here an open question, which should deserve further treatments according to the procedures defined in refs 57 and 58. Moving to benzoperiacenes, that is, nanoribbons of similar lengths but with a same number n of rings in each rank (see Figure 1), leads to the results of Table 3. In contrast to the previous series, for n = 4 (benzoperitetracene or dibenzoovalene), there is here no symmetry breaking, in agreement with our topological derivation. The symmetry breaking appears for n = 5 ( = 0.90), as expected from our model. For n = 7 a second symmetry breaking does not take place yet after left−right localization. Rhombuses. We have considered the series of rhombuses of increasing sizes, labeled according to the number n of rings in the width (short diagonal). As can be seen in Table 4, the model correctly predicts the lack of symmetry breaking for n = 2 and n = 3. A symmetry breaking appears in UDFT calculations for n = 4 (⟨S2⟩ = 0.79), in agreement with the model, the unpaired electrons tending to localize on the apex regions.87 For n = 5, the HOMO-1/LUMO+1 does not exhibit an intrinsic symmetry breaking, but after left−right localization of the two highest occupied orbitals and their mirror antibonding counterparts, the system should exhibit four unpaired electrons. This is consistent with the large value of ⟨S2⟩ observed in the UDFT calculation (⟨S2⟩ = 1.50), together with the calculated spin densities. They actually reveal two α electrons on the lateral atoms of one triangular subset and two β electrons on the lateral atoms of the other triangular subset, in consistency with scheme 3. Going to larger sizes might

seen as coronene substituted by two trimethylene-methane units (TMM), each bearing a spin one. For n = 4, the model again predicts a second symmetry breaking, while the corresponding calculated ⟨S2⟩ value now amounts to 2.06. This system can also be viewed as a tetrabenzo-ovalene carrying two TMM units, as discussed in ref 42. Linear Polyacenes. Linear polyacenes are usually seen as closed-shell systems. Nevertheless it is expected that increasing the length of the chain reduces the gap at the Fermi level and that this may induce specific effects. Performing UDFT calculations with B3LYP exchange potentials, Bendikov et al. found that a symmetry breaking takes place in hexacene and larger analogues.47 This finding, confirmed by several similar calculations,48−50 led to attribute diradical character to long polyacenes. It might however be considered as questionable.89 A significantin our mind decisiveanswer has been brought by Hachmann et al., who succeeded to perform a full valence π complete active space self-consistent field calculation (CASSCF), using the density matrix renormalization group technique.90 They calculated occupation numbers of the MOs (roots of one-particle density matrix) of these highly correlated functions. The evolution of the occupation numbers when increasing the chain length was plotted. The natural MOs where absolutely similar in shape to the CASSCF MOs, and they actually look extremely similar to the Hückel MOs. The occupation number of the nth MO, analogous to the HOMO, decreases when increasing the number of fused rings, while the occupation number of the next natural MO, similar to the LUMO, increases. Moreover, the same phenomenon takes place, although only for longer polyacenes, for the natural MOs resembling to the HOMO-1 and LUMO+1, respectively. It may be interesting to compare these high-quality results to those of our simple model. As shown in Table 6, the short polyenes up to heptacene (n = 7) should not present a symmetry breaking according to our model. Octacene (n = 8) should exhibit a symmetry breaking if U > 1.55|t|, the instability condition being almost satisfied for the ratio adopted in the present work (U = 1.5|t|). Actually, this molecule is the first polyacene exhibiting a symmetry breaking with our exchange-correlation potential and basis set (B3LYP/ 6-311G**). As spin densities on external rings are weak, these rings may be considered as essentially aromatic. In the central part of the molecule, two unpaired electrons are delocalized in undecapentaenyl tracks along the upper and lower chains (Scheme 4). It is interesting to notice that heptacene itself, and

Scheme 3. Best Clar’s Configuration in a Tetraradical Rhombus of Width 5

introduce larger numbers of unpaired electrons, but the modeling would require to proceed to more subtle localization of the highest occupied and lowest empty orbitals. Table 4 reports the evolution of the energy difference between the two singly occupied MOs of the ms = 1 UDFT solution. It decreases regularly and does not reveal by itself the appearance of several symmetry breakings in large enough architectures. Spiked Rhombuses. One may also consider the spiked rhombuses, or p-dimethylene rhombuses, in which extracyclic CH2 groups are added on apex atoms (see Figure 1, bottom).42 The results are reported in Table 5. For n = 2 a symmetry breaking takes place (⟨S2⟩ = 0.80) in agreement with the model. The molecule may be seen as a dimethylene pyrene (in agreement with the localization of the unpaired electrons, given by the pseudo spin densities of the UDFT calculation).42 After rotation of the MOs, our model predicts that a second symmetry breaking occurs for n = 3; this would be in line with the rather large value of ⟨S2⟩ (1.30), and the molecule can be

Scheme 4. Suggested Diradical Depiction of Octacene

all longer polyacenes, present a spontaneous dimerization,91 the covalent CC bonds between the monomers connecting the central atoms, consistently with their diradical picture. Increasing polyacene length does not seem to exhibit intrinsic instability with respect to (HOMO-1)/(LUMO+1) mixing. Yet, for lengths of 14 or 15 rings, the left−right localization of the two highest occupied and lowest unoccupied orbitals enables one to identify a double symmetry breaking and consequently to see the molecule as bearing four unpaired electrons. This is consistent with both the calculated value of J


Article

The Journal of Physical Chemistry A Scheme 5. Suggested Tetraradical Depiction of Decapentacene

⟨S2⟩ (larger than two for n = 14 and 15, see Table 6) and with the picture previously suggested for octacene. In longer polyacenes, in view of the spin densities one may say that two unpaired electrons of the same spin are spread along the upper edge, and two electrons of the other spin are spread along the lower edge. In a qualitative picture, it is tempting to see each of these pairs as two polyenyl monoradicals coupled ferromagnetically via a meta benzene ring. In this way, decapentacene (n = 15) can be viewed as a mixing and coupling of two octacene residues, as schematized in Scheme 5. This is actually supported by the calculated spin densities, exhibiting a clear depletion on median carbon atoms. The energy difference between the two singly occupied MOs of the ms = 1 UDFT calculation continuously decreases along the series up to n = 12 (see Table 6). But one may notice that for n > 12 this energy difference starts increasing (see Table 6). This is consistent with the fact that if the system is of 4unpaired electron character, the ms = 1 solution necessarily introduces a spin frustration. The relevant reference calculation concerns a ms = 2 calculation and the gap between the energies of the highest and lowest singly occupied MOs, which exhibits a regular decrease (see Table 6). The system exhibits a double symmetry breaking, as predicted by our model and in agreement with the occupation numbers provided by the full π CASSCF calculation of Hachmann et al.90 3.3. Comparison with UHF Instabilities. As previously mentioned the |U/t| ratio being much larger in HF calculations than for B3LYP-DFT, one may expect the occurrence of symmetry breaking for smaller members of the preceding series. Table 7 reports some significant results. The para-quinodi-

lowering (1.7 kcal/mol) but introduces a spin-density alternation along the ring, with rather large atomic spin densities (±0.61). In view of this result one may expect that UHF solutions will systematically appear in polycyclic systems. In the rhombus series pyrene itself is HF-unstable, with ⟨S2⟩ = 2.07, suggesting 4-unpaired electrons, but the spin density looks like a spin density wave, with spin-up densities on atoms of one color and spin-down densities on the other-color atoms which belong to the 0.79−0.85 interval. When increasing the size of the rhombus the ⟨S2⟩ value increases regularly, the mean value of the atomic spin densities being 0.90 for the [3]-rhombus and 0.93 for the [4]-rhombus. A qualitative view in terms of a few singly occupied MOs is difficult to suggest. In the series of polyacenes, naphthalene is HF-unstable with a value of ⟨S2⟩ (1.16) suggesting a diradical. However, by simply mixing the HOMO and the LUMO, one would not produce spin density on the atoms of the central bond, while the UHF spin densities again suggest a spin wave with ±0.77 spin density on all atoms. Actually the effective U/|t| ratio in HF is so large that one may think the electronic order of the π electron population from a magnetic viewpoint, that is, from a Heisenberg effective Hamiltonian. The corresponding model Hamiltonian only considers the neutral VB determinants and simply plays with the spin distributions, the effect of the ionic VB components goes through effective exchange coupling between adjacent atoms. This model Hamiltonian has proved to be a quantitative tool for the study of the low-energy covalent states of conjugated hydrocarbons.92 In this approach the leading valence-bond determinants are the so-called Néel determinants, presenting a full spin alternation between the atoms of both colors (hence on each chemical bond). The other spin distributions, which introduce some spin frustrations between adjacent atoms, have lower coefficients. The strong spin densities observed in the UHF calculations indicate a large collective effect, a huge spin polarization. In the large polycyclic systems the spin-symmetry breaking cannot be viewed as essentially governed by the mixing between the MOs close to the Fermi level, as was relevant in DFT calculations. The symmetry breaking is so strong in UHF calculation that starting from this solution to perform UMPn perturbative expansion of the wave function does not restore a proper S2 eigenfunction before very high orders.93 Special treatments are required to avoid this problem.94

Table 7. Mean Absolute Values of Alternating Spin Densities on Carbon Atoms in the BS ms = 0 Solutions Obtained from UHF/6-311G** Treatments atom spin densities ⟨S ⟩

mean

st. dev.

range

0.45 1.38 1.16 2.71 2.07 4.24 6.89

0.61 0.90 0.77 0.86 0.83 0.90 0.93

0.12 0.02 0.05 0.02 0.05 0.06

0.21 0.03 0.14 0.06 0.14 0.23

2

benzene p-quino-dimethane naphthalene perylene rhombus [2] (pyrene) rhombus [3] rhombus [4]

4. CONCLUSION The present work proposes a criterion which makes possible the identification of the number of occupied/virtual orbital pairs around the Fermi level, the mixing of which leads to intrinsic symmetry breaking(s). Accepting the idea that twice this number gives the number of open shells, this simple analysis offers a response to the question of the number of unpaired electrons in di- or polyradical molecules of the singlet ground state. The definition of the diradical character of a singlet state is somewhat arbitrary. In a rigorous way of thinking there is no qualitative difference between a closed-shell molecule and a

methane, as the first member of p-dimethylene polyphenylenes and of spiked rhombuses, already exhibits a strong symmetry breaking, with a large ⟨S2⟩ value (1.38) and a 30-kcal/mol stabilization with respect to the closed-shell solution. The spin densities range from 1.03 on the extra cyclic carbon to ±0.82 in the ring. In the periacene series, while the symmetry breaking takes place for the 4-length member, the perylene (length 2) is already unstable in HF, with a value of ⟨S2⟩ = 2.71, suggesting the existence of 4- or 6-unpaired electrons. But the mean spin density is 0.86, with a spread of 0.14. A surprising feature is the existence of a symmetry breaking in the benzene ring, which does not bring a large energy K


Article

The Journal of Physical Chemistry A

a tetra-radical character for this molecule. Of course, once the intrinsic symmetry breaking took place, localizing to some extent the two electrons involved in that symmetry breaking in different regions of the molecule, one may gain an additional energy by spin polarizing the other electrons, which fit the broken-symmetry exchange field. Predictions of this rather crude model have been validated by confronting them to the results of UDFT calculations. These calculations have been performed with the B3LYP exchange-correlation potential and if one assumes a ratio U = 1.5|t|, the predictions issued from the Hückel−Hubbard model Hamiltonian are verified on six series of polycyclic hydrocarbons. Of course the here-presented model is restricted to (i) alternant conjugated hydrocarbons (since it exploits the mirror theorem) and (ii) presenting small bond lengths variations, hence the choice of polycyclic systems. One might sophisticate it by using a geometry-dependent Hubbard Hamiltonian,63 but one would miss the simplicity of the predictive model. A second application of our analysis concerns a posteriori correction to the broken-symmetry mean-field solution, in order to reconstruct the CI matrix and get rid of the spincontamination problem. Knowing the energy difference between the closed-shell solution and the broken-symmetry solution, knowing the mean value of S2 for this solution, the quantities t and U can be evaluated from eqs 3 and 4, and then one is able to apply eq 6 to get a spin noncontaminated energy. By way of a last subtlety, let us recall that the same analysis might be applied to molecules presenting a high-spin ground state, for example, a triplet ground state with two singly occupied orbitals, with closed-shell levels possibly exhibiting intrinsic spin-symmetry breaking. In that case, its separation with spin polarization would be less obvious, though possible. As a general remark we would point out the intellectual benefit of the simple model Hamiltonians proposed 80 years ago for the treatment of electron systems in conjugated hydrocarbons and of general theorems established in the further decades by some quantum chemists as Longuett-Higgins. These theorems no longer appear in Quantum Chemistry textbooks, but they have predictive potentiality and should remain intellectual guides in the study of currently fashionable structures such as graphene flakes or nanotubes.

singlet di- or multiradical. Referring to the exact wave function, and looking at occupation numbers of natural orbitals, one may decide that the number of unpaired electrons is equal to the number of natural orbitals the occupation numbers of which are between 1.5 and 0.5, but these frontiers remain somewhat arbitrary. Being qualitative, the symmetry-breaking phenomenon presents the advantage of defining a frontier and the spin instability may be considered as a criterion for the existence of unpaired electrons. The mean value of the S2 operator brings an additional information. In principle, this mean value should be between 0 and 1 for two open shells and between 1 and 2 for four open shells. However, the spin polarization effect affects the ⟨S2⟩ value. A rigorous procedure follows the technique proposed by Coulaud et al.,57,58 which distinguishes spin instability and core spin polarization. If the ⟨S2⟩ value of a BS-UDFT solution is significantly larger than 1, this may come either from important spin polarization of the core by the two unpaired electrons or from the occurrence of a second symmetry breaking, resulting in the appearance of four unpaired electrons. Imposing a closed-shell character to n − 2 orbitals, one may check whether a second symmetry breaking has taken place and is responsible for this large observed value of ⟨S2⟩. This exploration must be performed, accepting eventually a simultaneous space-symmetry breaking of the MOs, which localizes them in disjoint regions of space, as illustrated in the numerical studies. The present work is devoted to conjugated polycyclic hydrocarbons, the ground state of which is a singlet, in agreement with Ochinnikov’s rule. By accepting Lewis-type pairing of electrons on bonds, they can be considered as having a closed-shell ground state. Nevertheless the mean-field treatments of these systems frequently show instability of the closed-shell solution, with a lower-energy solution appearing, in which α- and β-spin electrons occupy different space orbitals. The so-called spin-symmetry breaking, indicating a strong correlation effect and involving essentially the MOs around the Fermi levels, may take place in HF and KS-DFT calculations. In the HF case, however, the symmetry breaking is very strong, and it cannot be reduced to a mixing between frontier orbitals and behaves as a collective phenomenon. The present work shows that one may predict at a low cost whether a given molecular architecture will be subject to such symmetry breaking in KS-DFT calculations. The prediction is based on essentially topological model Hamiltonians, namely, (i) the tight-binding Hückel Hamiltonian, with uniform on-site energies and unique hopping integral t between adjacent carbon atoms, and (ii) the extremely simplified Hubbard bielectronic Hamiltonian, reduced to on-site repulsion U. Provided that a |t|/U ratio has been determined as corresponding to the chosen exchange-correlation used in KS-DFT calculations, one may predict from the Hückel solution whether the mean-field treatment of any conjugated hydrocarbon will exhibit a spin-symmetry breaking. The present analysis has shown that the condition for symmetry-breaking occurrence may be formulated very simply in terms of the energy of an occupied MO and of the repulsion between the two electrons occupying this closed-shell MO. As soon as their ratio is smaller than 1, the symmetry breaking takes place, and the molecule may then be considered as a singlet diradical. The first orbital candidate for symmetry breaking is of course the HOMO/LUMO couple, but in some systems the HOMO-1/LUMO+1 couple of orbitals may undergo intrinsic symmetry breaking as well, then suggesting

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b07597. CI transcription of spin polarization in symmetryadapted formalisms (PDF)

■

AUTHOR INFORMATION

ORCID

Jean-Paul Malrieu: 0000-0003-0868-8391 Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Overhauser, A. W. Structure of nuclear matter. Phys. Rev. Lett. 1960, 4, 415−418. (2) Thouless, D. J. In Quantum Mechanics of Many-Body Systems; Academic Press: New York, 1961. (3) Adams, W. H. Stability of Hartree-Fock states. Phys. Rev. 1962, 127, 1650−1658. L


Article

The Journal of Physical Chemistry A (4) Löwdin, P. O. Discussion on the Hartree-Fock approximation. Rev. Mod. Phys. 1963, 35, 496−501. (5) Koutecky, J. Unrestricted Hartree-Fock solutions for closed-shell molecules. J. Chem. Phys. 1967, 46, 2443−2445. (6) Lefebvre, R.; Smeyers, Y. Extended Hartree-Fock calculations for the helium ground state. Int. J. Quantum Chem. 1967, 1, 403−419. (7) Cizek, J.; Paldus, J. Stability conditions for the solutions of the Hartree-Fock equations for atomic and molecular systems. Application to the p-electron model of cyclic polyenes. J. Chem. Phys. 1967, 47, 3976−3985. (8) Paldus, J.; Cizek, J. Stability conditions for the solutions of the Hartree-Fock equations for atomic and molecular systems.V. The nonanalytic behavior of the broken-symmetry solutions at the branching point. Phys. Rev. A: At., Mol., Opt. Phys. 1971, 3, 525−527. (9) Noodleman, L. Valence bond description of antiferromagnetic coupling in transition metal dimers. J. Chem. Phys. 1981, 74, 5737− 5743. (10) Lepetit, M. B.; Malrieu, J.-P.; Trinquier, G. Valence bond formulation of Hartree-Fock instability conditions for simple and multiple bonds. Chem. Phys. 1989, 130, 229−239. (11) Yamaguchi, K.; Fukui, H.; Fueno, T. Molecular orbital (MO) theory for magnetically interacting organic compounds. Ab-initio MO calculations of the effective exchange integrals for cyclophane-type dimers. Chem. Lett. 1986, 15, 625−628. (12) Yamaguchi, K.; Takahara, Y.; Fueno, T.; Houk, K. N. Extended Hartree-Fock (EHF) theory of chemical reactions. Theor. Chim. Acta 1988, 73, 337−364. (13) Yamanaka, S.; Okumura, M.; Nakano, M.; Yamaguchi, K. EHF theory of chemical reactions. Part 4. UNO CASSCF, UNO CASPT2 and R(U)HF coupled-cluster (CC) wavefunctions. J. Mol. Struct. 1994, 310, 205−218. (14) Caballol, R.; Castell, O.; Illas, F.; de P. R. Moreira, I.; Malrieu, J.P. Remarks on the proper use of the broken symmetry approach to magnetic coupling. J. Phys. Chem. A 1997, 101, 7860−7866. (15) Ferré, N.; Guihéry, N.; Malrieu, J.-P. Spin decontamination of broken-symmetry density functional theory calculations: deeper insight and new formulations. Phys. Chem. Chem. Phys. 2015, 17, 14375−14382. (16) Bulik, I. W.; Henderson, T. M.; Scuseria, G. Can singlereference coupled cluster theory describe static correlation ? J. Chem. Theory Comput. 2015, 11, 3171−3179. (17) Colvin, M. E.; Janssen, C. L.; Seidl, E. T.; Nielsen, I. B. M.; Melius, C. F. Energies, resonance and UHF instabilities in polycyclic aromatic hydrocarbons and linear polyenes. Chem. Phys. Lett. 1998, 287, 537−541. (18) Bozkaya, U. Analytic energy gradients and spin multiplicities for orbital-optimized second-order perturbation theory with density-fitting approximation: an efficient implementation. J. Chem. Theory Comput. 2014, 10, 4389−4399. (19) Lepetit, M. B.; Malrieu, J.-P.; Pelissier, M. Multiplicity of symmetry-broken Hartree-Fock solutions in multiple bonds and atomic clusters: an asymptotic view. Phys. Rev. A: At., Mol., Opt. Phys. 1989, 39, 981−991. (20) Cohen, R. D.; Sherrill, C. D. The performance of density functional theory for equilibrium molecular properties of symmetry breaking molecules. J. Chem. Phys. 2001, 114, 8257−8269. (21) Bozkaya, U. Orbital-optimized third-order Møller-Plesset perturbation theory and its spin- component and spin-opposite scaled variants: application to symmetry breaking problems. J. Chem. Phys. 2011, 135, 224103−224113. (22) Bozkaya, U.; Sherrill, C. D. Analytic energy gradients for the orbital-optimized second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2013, 138, 184103−184108. (23) McAdon, M. H.; Goddard, W. A., III Charge density waves, spin density waves, and Peierls distortions in one-dimensional metals. I. Hartree-Fock studies of Cu, Ag, Au, Li, and Na. J. Chem. Phys. 1988, 88, 277−302. (24) Lepetit, M. B.; Malrieu, J.-P.; Spiegelmann, F. Toward a magnetic description of metals in terms of interstitial molecular

orbitals: exploiting the multiplicity of symmetry-broken Hartree-Fock solutions on small alkali-metal clusters. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 8093−8106. (25) Ovchinnikov, A. A. Multiplicity of the ground state of large alternant organic molecules with conjugated bonds. Theor. Chim. Acta 1978, 47, 297−304. (26) Lieb, E. H. Two Theorems on the Hubbard Model. Phys. Rev. Lett. 1989, 62, 1201−1204. (27) Malrieu, J.-P.; Ferré, N.; Guihéry, N.Magnetic properties of conjugated hydrocarbons from topological Hamiltonians. In Challenges in Computational Chemical Physics, Vol. 22; Chauvin, R., Lepetit, C., Silvi, B., Alikani, E., Eds; Applications of Topological Methods in Molecular Chemistry; Springer: Heidelberg, 2016; pp 361−395. (28) See for instance Carbon-based magnetism: an overview of metalfree carbon-based compounds and materials; Makarova, T. L., Palacio, F., Eds.; Elsevier: Amsterdam, 2006. (29) Rajca, S.; Rajca, A.; Wongsriratanakul, J.; Butler, P.; Choi, S.-M. Organic spin clusters. A dendritic-macrocyclic poly(arylmethyl)polyradical with very high spin of S = 10 and its derivatives: synthesis, magnetic studies, and small-angle neutron scattering. J. Am. Chem. Soc. 2004, 126, 6972−6986. (30) Borden, W. T. Diradicals; Wiley-Interscience: New York, 1982. Du, P.; Borden, W. T. Ab initio calculations predict a singlet ground state for tetramethyleneethane. J. Am. Chem. Soc. 1987, 109, 930−931. (31) Nicolaides, A.; Borden, W. T. CI calculations on didehydrobenzenes predict heats of formation for the meta and para isomers that are substantially higher than previous experimental values. J. Am. Chem. Soc. 1993, 115, 11951−11957. (32) Borden, W. T.; Iwamura, H.; Berson, J. A. Violations of Hund’s rule in non-Kekulé hydrocarbons: theoretical prediction and experimental verification. Acc. Chem. Res. 1994, 27, 109−116. (33) Hrovat, D. A.; Borden, W. T. MCSCF and CASPT2N calculations on the excited states of 1,2,4,5-tetramethylenebenzene. The UV-vis spectrum observed belongs to the singlet state of the diradical. J. Am. Chem. Soc. 1994, 116, 6327−6331. (34) Borden, W. T.; Davidson, E. R.; Feller, D. RHF and twoconfiguration SCF calculations are inappropriate for conjugated diradicals. Tetrahedron 1982, 38, 737−739. Borden, W. T.; Davidson, E. R. The importance of including dynamic electron correlation in ab Initio calculations. Acc. Chem. Res. 1996, 29, 67−75. (35) Fang, S.; Lee, M.-S.; Hrovat, D. A.; Borden, W. T. Ab-initio calculations show why m-phenylene is not always a ferromagnetic coupler. J. Am. Chem. Soc. 1995, 117, 6727−6731. (36) Bally, T.; Borden, W. T. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; Wiley-VCH: New York, 1999, Vol. 13; pp 1−97. (37) Lineberger, W. C.; Borden, W. T. The synergy between qualitative theory, quantitative calculations, and direct experiments in understanding, calculating, and measuring the energy differences between the lowest singlet and triplet states of organic diradicals. Phys. Chem. Chem. Phys. 2011, 13, 11792−11813. (38) Abe, M. Diradicals. Chem. Rev. 2013, 113, 7011−7088. (39) Filatov, M.; Shaik, S. Tetramethyleneethane (TME) diradical: experiment and density functional theory reach an agreement. J. Phys. Chem. A 1999, 103, 8885−8889. (40) Rodríguez, E.; Reguero, M.; Caballol, R. The controversial ground state of tetramethyleneethane. An ab initio CI study. J. Phys. Chem. A 2000, 104, 6253−6258. (41) Pozun, Z. D.; Su, X.; Jordan, K. D. Establishing the ground state of the disjoint diradical tetramethyleneethane with quantum Monte Carlo. J. Am. Chem. Soc. 2013, 135, 13862−13869. (42) Trinquier, G.; Malrieu, J.-P. Kekulé versus Lewis: when aromaticity prevents electron pairing and imposes polyradical character. Chem. - Eur. J. 2015, 21, 814−828. ̈ (43) Thiele, J.; Balhorn, H. Ueber einen chinoiden kohlenwasserstoff. Ber. Dtsch. Chem. Ges. 1904, 37, 1463−1470. (44) Tschitschibabin, A. E. Ü ber einige phenylierte derivate des p, pditolyls. Ber. Dtsch. Chem. Ges. 1907, 40, 1810−1819. Montgomery, L. K.; Huffman, J. C.; Jurczak, E. A.; Grendze, M. P. The molecular M


Article

The Journal of Physical Chemistry A structures of Thiele’s and Chichibabin’s hydrocarbons. J. Am. Chem. Soc. 1986, 108, 6004−6011. (45) Muller, E.; Pfanz, H. Ü ber biradikaloide terphenylderivate. Ber. Dtsch. Chem. Ges. B 1941, 74, 1051−1074. (46) Shishlov, N. M.; Asfandiarov, N. L. Russ. Chem. Bull. 2000, 49, 1676−1681. (47) Bendikov, M.; Duong, H. M.; Starkey, K.; Houk, K. N.; Carter, E. A.; Wudl, F. Oligoacenes: theoretical prediction of open-shell singlet diradical ground states. J. Am. Chem. Soc. 2004, 126, 7416− 7417. (48) Dos Santos, M. C. Electronic properties of acenes: oligomer to polymer structure. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 045426−9. (49) Jiang, D.; Dai, S. Electronic ground state of higher acenes. J. Phys. Chem. A 2008, 112, 332−335. (50) Rayne, S.; Forest, K. A comparison of density functional theory (DFT) methods for estimating the singlet−triplet (S0−T1) excitation energies of benzene and polyacenes. Comput. Theor. Chem. 2011, 976, 105−112. (51) Torres, A. E.; Guadarrama, P.; Fomine, S. Multiconfigurational character of the ground states of polycyclic aromatic hydrocarbons. A systematic study. J. Mol. Model. 2014, 20, 2208−2216. (52) Hückel, E. Zur quantentheorie der doppelbindung. Eur. Phys. J. A 1930, 60, 423−456. Hückel, E. Quantentheoretische beiträge zum benzolproblem. Eur. Phys. J. A 1931, 70, 204−286. Hückel, E. Zur quantentheorie der doppelbindung und ihres stereochemischen verhaltens. Z. Elektrochem. Angew. Phys. Chem. 1930, 36, 641−645. (53) Hubbard, J. Electron correlations in narrow energy bands. Proc. R. Soc. London, Ser. A 1963, 276, 238−257. (54) Malrieu, J.-P.; Trinquier, G. Proper use of broken-symmetry calculations in antiferromagnetic polyradicals. J. Chem. Phys. 2016, 144, 211101−3. (55) Longuet-Higgins, H. C. Some studies in molecular orbital theory. I. Resonance structures and molecular orbitals in unsaturated hydrocarbons. J. Chem. Phys. 1950, 18, 265−274. (56) Hall, G. G. The bond orders of alternant hydrocarbon molecules. Proc. R. Soc. London, Ser. A 1955, 229, 251−259. (57) Coulaud, E.; Guihéry, N.; Malrieu, J.-P.; Hagebaum-Reignier, D.; Siri, D.; Ferré, N. Analysis of the physical contributions to magnetic couplings in broken-symmetry density functional theory approach. J. Chem. Phys. 2012, 137, 114106−8. (58) Coulaud, E.; Malrieu, J.-P.; Guihéry, N.; Ferré, N. Additive decomposition of the physical components of the magnetic coupling from broken symmetry density functional theory calculations. J. Chem. Theory Comput. 2013, 9, 3429−3436. (59) Stück, D.; Baker, T. A.; Zimmerman, P.; Kurlancheek, W.; Head-Gordon, M. On the nature of electron correlation in C60. J. Chem. Phys. 2011, 135, 194306−5. (60) Edmiston, C.; Ruedenberg, K. Localized atomic and molecular orbitals. Rev. Mod. Phys. 1963, 35, 457−465. (61) Gaussian 09, Revision D.01; Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., Scalmani, G., Barone, V., Mennucci, B., Petersson, G. A., et al. Gaussian, Inc.: Wallingford, CT, 2009. (62) Full geometry optimizations were carried out up to energy gradients lower than 10−5 au. Cartesian coordinates of optimized geometries are available upon request. (63) Garcia, V. M.; Caballol, R.; Malrieu, J.-P. Theoretical study of the ethylene electronic spectrum and extraction of an r-dependent Hubbard Hamiltonian. Chem. Phys. Lett. 1996, 261, 98−104. (64) Calzado, C. J.; Cabrero, J.; Malrieu, J.-P.; Caballol, R. Analysis of the magnetic coupling in binuclear complexes. II. Derivation of valence effective Hamiltonians from ab initio CI and DFT calculations. J. Chem. Phys. 2002, 116, 3985−4000. (65) Schmalz, T. G.; Serrano-Andres, L.; Sauri, V.; Merchan, M.; Oliva, J. M. A distance-dependent parameterization of the extended Hubbard model for conjugated and aromatic hydrocarbons derived from stretched ethene. J. Chem. Phys. 2011, 135, 194103−10.

(66) In theoretical PAH databases, the ir spectrum for this neutral C36H16 molecule is usually obtained from restricted DFT treatment of the closed-shell configuration (cf. for instance: Malloci, G.; Joblin, C.; Mulas, G. On-line database of the spectral properties of polycyclic aromatic hydrocarbons. Chem. Phys. 2007, 332, 353−359. (67) Clar, E.Polycyclic Hydrocarbons; Academic Press: London, 1964; Vol. 2. Clar, E.The Aromatic Sextet; Wiley: New York, 1972. (68) Portella, G.; Poater, J.; Sola, M. Assessment of Clar’s aromatic πsextet rule by means of PDI, NICS and HOMA indicators of local aromaticity. J. Phys. Org. Chem. 2005, 18, 785−791. (69) Maksić, Z. B.; Barić, D.; Müller, T. Clar’s sextet rule is a consequence of the σ-electron framework. J. Phys. Chem. A 2006, 110, 10135−10147. (70) Balaban, A. T.; Klein, D. J. Claromatic carbon nanostructures. J. Phys. Chem. C 2009, 113, 19123−19133. (71) Wassmann, T.; Seitsonen, A. P.; Saitta, A. M.; Lazzeri, M.; Mauri, F. Clar’s theory, π-electron distribution, and geometry of graphene nanoribbons. J. Am. Chem. Soc. 2010, 132, 3440−3451. (72) Fujii, S.; Enoki, T. Clar’s aromatic sextet and π-electron distribution in nanographene. Angew. Chem., Int. Ed. 2012, 51, 7236− 7241. (73) Sola, M. Forty years of Clar’s aromatic π-sextet rule. Front. Chem. 2013, 1, 1−8. (74) Gutman, I.; Radenkovic, S.; Antic, M.; Djurdjevic, J. A test of Clar aromatic sextet theory. J. Serb. Chem. Soc. 2013, 78, 1539−1546. (75) For the use of Clar’s sextets in contexts closer to the present work, see Das, S.; Wu, J. in Polycyclic arenes and heteroarenes: synthesis, properties, and applications, 1st ed.; Miao, Q., Ed.; Wiley-VCH Verlag GmbH & Co, 2016; Chapter 1, p 3; “Open-shell benzenoid polycyclic hydrocarbons”. (76) Jiang, D.-E; Dai, S. Circumacenes versus periacenes: HOMO− LUMO gap and transition from nonmagnetic to magnetic ground state with size. Chem. Phys. Lett. 2008, 466, 72−75. (77) Moscardó, F.; San-Fabián, E. On the existence of a spinpolarized state in the n-periacene molecules. Chem. Phys. Lett. 2009, 480, 26−30. (78) Gutman, I.; Durdevic, J.; Radenkovic, S.; Matovic, Z. Anomalous cyclic conjugation in the perylene/bisanthrene homologous series. Monatsh. Chem. 2012, 143, 1649−1653. (79) Sun, Z.; Ye, Q.; Chi, C.; Wu, J. Low band gap polycyclic hydrocarbons: from closed-shell near infrared dyes and semiconductors to open-shell radicals. Chem. Soc. Rev. 2012, 41, 7857− 7889. (80) Sun, Z.; Wu, J. Open-shell polycyclic aromatic hydrocarbons. J. Mater. Chem. 2012, 22, 4151−4160. (81) Plasser, F.; Pasalic, H.; Gerzabek, M. H.; Libisch, F.; Reiter, R.; Burgdörfer, J.; Müller, T.; Shepard, R.; Lischka, H. The multiradical character of one- and two-dimensional graphene nanoribbons. Angew. Chem., Int. Ed. 2013, 52, 2581−2584. (82) Martin-Martinez, F. J.; Fias, S.; Van Lier, G.; De Proft, F.; Geerlings, P. Tuning aromaticity patterns and electronic properties of armchair graphene nanoribbons with chemical edge functionalisation. Phys. Chem. Chem. Phys. 2013, 15, 12637−12647. (83) Dias, J. R. Valence-bond determination of diradical character of polycyclic aromatic hydrocarbons: from acenes to rectangular benzenoids. J. Phys. Chem. A 2013, 117, 4716−4725. (84) Ye, Q.; Chi, C. Recent highlights and perspectives on acene based molecules and materials. Chem. Mater. 2014, 26, 4046−4056. (85) Sun, Z.; Wu, J. Closed-shell and open-shell 2D nanographenes. Top. Curr. Chem. 2012, 349, 197−248. (86) Sun, Z.; Zeng, Z.; Wu, J. Zethrenes, extended p-quinodimethanes, and periacenes with a singlet biradical ground state. Acc. Chem. Res. 2014, 47, 2582−2591. (87) Yeh, C.-N.; Chai, J.-D. Role of Kekulé and non-Kekulé structures in the radical character of alternant polycyclic aromatic hydrocarbons: A TAO-DFT study. Sci. Rep. 2016, 6, 30562. (88) Torres, A. E.; Flores, R.; Fomine, S. A comparative study of one and two dimensional π-conjugated systems. Synth. Met. 2016, 213, 78−87. N


Article

The Journal of Physical Chemistry A (89) Huzak, M.; Deleuze, M. S.; Hajgató, B. Half-metallicity and spincontamination of the electronic ground state of graphene nanoribbons and related systems: an impossible compromise ? J. Chem. Phys. 2011, 135, 104704−18. (90) Hachmann, J.; Dorando, J. J.; Avilés, M.; Chan, G. K.-L. The radical character of the acenes: a density matrix renormalization group study. J. Chem. Phys. 2007, 127, 134309−19. (91) Mondal, R.; Tonshoff, C.; Khon, D.; Neckers, D. C.; Bettinger, H. F. Synthesis, stability and photochemistry of pentacene, hexacene and heptacene: a matrix isolation study. J. Am. Chem. Soc. 2009, 131, 14281−14289. (92) Said, M.; Maynau, D.; Malrieu, J.-P.; Garcia-Bach, M. A. A nonempirical Heisenberg Hamiltonian for the study of conjugated hydrocarbons. Ground state conformational studies. J. Am. Chem. Soc. 1984, 106, 571−579. (93) Lepetit, M. B.; Pelissier, M.; Malrieu, J.-P. Origins of the poor convergence of Many-Body perturbation theory from Unrestricted Hartree-Fock zero-order descriptions. J. Chem. Phys. 1988, 89, 998− 1008. (94) Bozkaya, U. Analytic energy gradients and spin multiplicities for orbital-optimized second-order perturbation theory with density-fitting approximation: an efficient implementation. J. Chem. Theory Comput. 2014, 10, 4389−4399.

O


Can a Topological Approach Predict Spin-Symmetry Breaking in

Recommend Documents