A Recipe for Geometry Optimization of Diradicalar Singlet States from

Article pubs.acs.org/JPCA

A Recipe for Geometry Optimization of Diradicalar Singlet States from Broken-Symmetry Calculations Jean-Paul Malrieu* and Georges Trinquier Laboratoire de Chimie et Physique Quantiques (CNRS-UMR5626), IRSAMC, Université Paul-Sabatier (Toulouse III), 31062 Toulouse Cedex, France ABSTRACT: The equilibrium geometries of the singlet and triplet states of diradicals may be somewhat different, which may have an influence on their magnetic properties. The single-determinantal methods, such as Hartree−Fock or Kohn−Sham density functional theory, in general rely on broken-symmetry solutions to approach the singlet-state energy and geometry. An approximate spin decontamination is rather easy for the energy of this state but is rarely performed for its geometry optimization. We suggest simple procedures to estimate the optimized geometry and energy of a spin-decontaminated singlet, the accuracies of which are tested on a few organic diradicals. This technique can be generalized to interactions between higher-spin units or to multispin systems.

1. INTRODUCTION In its simple Kohn−Sham version, the density functional theory (DFT) is extremely convenient and popular. Of course, the exchange correlation potential is not known, a very large number of variants are practiced, most of them introducing some parameters, but a reasonable accuracy is attainable for nonexotic closed-shell molecules in their ground state, especially regarding the equilibrium geometries and vibrational frequencies. The use of DFT is also very popular for magnetic systems which are of open-shell character, with a few singly occupied orbitals.1 The state of highest spin multiplicity and highest Ms value may be approached by a single determinant, with 2 Ms electrons in 2 Ms singly occupied orbitals with parallel spins, and the single-determinantal approaches are in principle applicable to estimate its energy. One may use either restricted or unrestricted formalisms, the former one keeping a closed-shell character to the doubly occupied orbitals, the second one introducing some spin polarization, which means that slightly different orbitals are used for electrons of α and β spins. In this case, the determinant is not a pure spin state, i.e., an eigenfunction of the S2 operator, but the deviation of its mean value from its expected value is weak. If geometry optimization of so-described states of highest spin multiplicity does not present any conceptual difficulty, a severe problem concerns the states of lower spin multiplicities, which cannot be represented as an open-shell single determinant. One may optimize an open-shell single determinant of low Ms value, with Ms α and β unpaired electrons tending to localize in different sites. This solution is a broken-symmetry (BS) determinant,2−5 © 2012 American Chemical Society

and the corresponding function is a mixture of spin eigenfunctions. From the energies of the upper multiplet and BS determinants, it is possible to approximate the value of the lowest S2 eigenvalue, according to one of the more-or-less approximate spin-decontamination techniques.6−8 Regarding the optimal geometry of the low-spin state, it cannot in principle be obtained by minimizing the energy of the BS determinant, as noticed in various works where principles of appropriate geometry optimization of the low-spin state are proposed.9 The present work first suggests a simple analytical procedure to estimate the optimal geometry and energy of the low-spin state from the only knowledge of geometries and energies of the single-determinant solutions of high and low Ms. Next, a more refined treatment will exploit the optimized geometries of the triplet and broken-symmetry determinants to define a collective coordinate along which a spin-decontaminated energy curve is calculatedhence a better evaluation of singlet-state energy minimum. Both treatments are applied on a series of organic diradicals presenting either a ferromagnetic or an antiferromagnetic coupling between the unpaired electrons. Last, subsequent generalizations to any high-spin units will be assessed. Received: April 20, 2012 Revised: July 5, 2012 Published: July 27, 2012 8226

dx.doi.org/10.1021/jp303825x | J. Phys. Chem. A 2012, 116, 8226−8237

The Journal of Physical Chemistry A

Article

2. SITUATION OF THE PROBLEM A. Energy Difference. Let us specify the problem for the elementary case of two unpaired electrons in two molecular orbitals a and b, which may be localized in different parts of the molecule. If these orbitals are orthogonal, the triplet-state single determinant will be written as

Es − E T =

This formula has the advantage of being applicable whatever the diradicalar character of the singlet state, but it may present a defect in the present context. If the Ms=0 solution tends to become a closed shell, i.e., a pure singlet, the singlet-to-triplet energy difference is no longer equal to EBS − ET as it should. Due to the spin polarization of the core orbitals, ⟨S2⟩T is here not equal to 2, and its deviation also depends on the spin polarizationa real and physically meaningful phenomenon that is not to be ignored. Alternatively, one may solve the two equations regarding the mean values of H and S2 on the BS solution, which is a linear combination of the sought singlet state and of the Ms=0 solution for the triplet

Φ+T = |Πii i ̅ ab| in a restricted formalism, and Φ′+T = |Πii i ̅ ′ab|

in the unrestricted one. The mean values of the S2 operator is equal to 2 for the first determinant, and close to 2 for the second one. The i index runs over all doubly occupied orbitals. In the second solution, the difference between the “core” orbitals i and i′ introduces the effect of the exchange field created by the two unpaired electrons of α spin, in a and b, on the α-spin electrons of the core. This is a spin-polarization phenomenon, responsible, for instance, of the appearance of a σ spin density in planar conjugated hydrocarbons where unpaired electrons occupy π orbitals. Starting from a restricted openshell Fock operator, the operator responsible for spin polarization of doubly occupied orbitals is (Ka+Kb)/2, where Ka and Kb are the exchange operators associated to the singly occupied orbitals a and b. Let us call ET the energy of Φ′+T

ET =

2(E BS − E T) ⟨S2⟩T − ⟨S2⟩BS

ΦBS = γ ΦS + δ ΦT(Ms = 0) ⟨ΦBS|S2|ΦBS⟩ = δ 2⟨ΦT|S2|ΦT⟩

(assuming that the mean value of S2 is the same for the Ms=0 and Ms=1 components of the triplet) ⟨ΦBS|H |ΦBS⟩ = γ 2⟨ΦS|H |ΦS⟩ + δ 2⟨ΦT|H |ΦT⟩

which leads to Es − E T =

⟨Φ′+T |H |Φ′+T ⟩

The corresponding proper singlet state should be written as ΦS = |Πiii(ab̅ + ba)|/√2, i.e., from a two-determinant wave ̅ ̅ function. One may of course optimize the energy of a single determinant ΦBS = |Πii′i′(a′b̅′)|,10 or more generally, in an unrestricted formalism, ΦBS = |Πii″i‴(a′b̅′)|. Let us call EBS its energy

⟨S2⟩T (E BS − E T) ⟨S2⟩T − ⟨S2⟩BS

With this formula, the gap behaves correctly when the Ms=0 solution tends to become a closed shell. Note that the simple relation Es − E T =

E BS = ⟨ΦBS|H |ΦBS⟩

2(E BS − E T) 2 − ⟨S2⟩BS

also behaves correctly in this limiting situation and may as well be used in this case. At this point, it is convenient to introduce a general spin-decontamination factor (SDF), which can be taken at 2, or defined either from ref 6 as

The essentially singly occupied orbitals a′ and b′ have no reason to remain orthogonal, their delocalization introduces some valence-bond (VB) ionic components which are responsible for the stabilization of the singlet state, according to the kinetic exchange (or superexchange) mechanism.11 The difference between the essentially closed-shell orbitals, i″ and i‴, introduces the spin polarization of the core, the electrons of which feel the presence of an α-spin electron near a and of a βspin one near b, through an exchange field (Ka − Kb)/2.12 If the overlap between a′ and b′ remains small and if spin polarization is weak enough, the BS determinant may be considered as a half-and-half mixture of the desired singlet state ΦS and of the Ms=0 component of the triplet state

λ′ =

2 ⟨S ⟩T − ⟨S2⟩BS 2

or more generally as λ=

⟨S2⟩T 2

⟨S ⟩T − ⟨S2⟩BS

(1)

The corrected singlet-to-triplet energy difference therefore writes as Es − E T = λ(E BS − E T)

ΦT0 = |Πiii(ab ̅ − ba ̅ )| / 2

B. Geometries. The geometry optimization of the triplet state provides an energy minimum Emin and an optimal T geometry Xmin , corresponding to a set of internal coordinates T (bond lengths, bond angles and dihedral angles). The energy of the BS solution calculated for this geometry,

The corresponding mean value of the S2 operator ⟨S2⟩BS = ⟨ΦBS|S2|ΦBS⟩

is then close to 1. Its energy should be halfway between the energy of the triplet state and that of the singlet state, so that the singlet-state energy may be written as

* = E BS(X Tmin) E BS leads to the vertical transition energy from the optimal triplet geometry (Figure 1)

Es = E T + 2(E BS − E T) ⇒ Es − E T = 2(E BS − E T)

A more sophisticated estimate of the singlet−triplet energy difference has been proposed6,7 by Yamaguchi as

* − E Tmin) ES* − E Tmin = λ(E BS 8227



Article

Figure 1. Schematic representation of potential energy along the collective geometrical coordinate Q. min = 1), λ = 2, therefore locating Qmin S around 2. For Q = QS , the optimized energy of the singlet state is

One may also determine the optimal geometry Xmin BS of the BS determinant and its energy Emin BS . The energy of the triplet for this geometry is

* − E Tmin) − ESmin = E Tmin + λ(E BS

min E T* = E T(XBS )

Of course δET=ET*−Emin and δEBS=EBS * −Emin T BS are positive. If ⟨S2⟩BS ≠ 0, the geometry Xmin is not the optimal geometry of BS the spin-decontaminated singlet state, but one may assume that min the Xmin T → XBS geometry change goes in the right direction and that the spin-decontaminated singlet geometry is lying somewhere along this collective geometrical change. Let us focus on this collective motion path, quantified by an advancement parameter Q:

which gives the adiabatic triplet-to-singlet energy gap * − E Tmin) − ESmin − E Tmin = λ(E BS

Again, this expression can take simplified forms. If λ = 2, it becomes * − E Tmin) − ESmin − E Tmin = 2(E BS

Assuming an harmonic behavior of the potential energy surfaces, the quantities δET and δEBS are sufficient to determine the curvatures of the curves ET(Q) and EBS(Q):

* − E Tmin) − 4δE BS ESmin − E Tmin = 2(E BS

If the spin-decontamination factor λ does not depend on Q, the singlet-state energy may be written as

(6)

* − E Tmin) − λ 2δE BS ESmin − E Tmin = λ(E BS

min ES(Q ) = λ(E BS + (Q − 1)2 δE BS)

All these derivations suppose that the decontamination factor λ does not change significantly with the geometry. If this parameter is different for the Ms=1 and Ms=0 optimized geometries one may either (i) try to use a linear dependence of λ on the geometry change,

(2)

Direct differentiation provides the optimized geometry of the singlet state λδE BS

(5)

suggesting a singlet-state relaxation from triplet geometry four times larger than that of the BS solution. Still assuming δEBS = δET, but now taking λ ≠ 2, one gets

min E BS(Q ) = E BS + (Q − 1)2 δE BS

λδE BS − (λ − 1)δE T

4(δE BS)2 2δE BS − δE T

If both triplet and BS potential curves exhibit a same curvature (δEBS = δET), it further reduces to

E T(Q ) = E Tmin + Q 2δE T

Q Smin =

(λδE BS)2 λδE BS − (λ − 1)δE T (4)

min X(Q ) = X Tmin + Q (XBS − X Tmin)

+ (1 − λ)(E Tmin + Q 2δE T)

(λδE BS)2 λδE BS − (λ − 1)δE T

λ(Q ) = λ 0 + λ1Q (λ0 and λ1 being the values of λ at geometries of triplet and BS solutions, respectively) or, preferably, (ii) recalculate the energies of the triplet and BS determinants along the Q coordinate. This numerical method consists in recalculating the

(3)

If the curvatures are identical (δEBS = δET), the above relation 2 2 2 simplifies to Qmin S = λ. For ideal values of ⟨S ⟩ (⟨S ⟩T = 2, ⟨S ⟩BS 8228



Article

spin-decontaminated energy for a series of values of Q between 0 and 2 (and beyond if necessary), λ being recalculated for each value of Q: Es(Q ) = E T(Q ) + λ(Q )(E BS(Q ) − E T(Q ))

(7)

This method consists of a constrained exploration of the potential-energy surface of the decontaminated singlet. Due to this constraint, it only leads to an upper bound of the overall energy minimum. This technique is somewhat more expensive, but it offers the possibility to check the validity of the predictions based on the preceding analytical formulas.

3. NUMERICAL TESTS In coordination-chemistry magnetic complexes, when going from triplet to singlet state of a diradical, geometry changes are expected to be extremely small, and are usually neglected. Actually in these systems the unpaired electrons are well localized on the transition metal ions, connected through rather rigid ligands. As geometry changes are supposed to be larger in organic diradicals, especially in conjugated hydrocarbons, where unpaired electrons are spread on the entire π system, it is therefore desirable to study this phenomenon on organic systems. Because one may suspect that equilibrium-geometry changes are larger when the magnetic coupling (or singlet-totriplet transition energy) is large, the selected systems should also present non-negligible singlet−triplet differences. Three examples have been selected. The first one is a typical ferromagnetic system, m-xylylene, 1, the ground state of

which is triplet, with a transition energy to the singlet state around 0.5 eV.13 As an antiferromagnetic coupling, we chose pbenzyne, 2, in which the two unpaired electrons belong to the σ system, and where the singlet-to-triplet transition energy is around 0.3 eV.14 Last, we selected an architecture constituted by two phenalene radicals coupled in an antiferromagnetic manner by an acetylene bridge, 3, a system already studied by Nakano et al. and possibly presenting dramatic properties.15 All the calculations were performed using the UDFT version of the

Figure 2. Energy curves along the collective geometrical coordinate Q for m-xylylene. 8229



Article

Table 1. Calculated QMin along the Collective Geometrical S Coordinate

Gaussian03 quantum chemistry package, with standard B3LYP hybrid functional and 6-311G** basis set.16 Optimized geometries were carried out up to energy gradients lower than 10−5 au. Optimized geometries correspond to C2v, D2h, and C2h symmetries in the three molecules, respectively. If the sensitivity of magnetic coupling to the percentage of exact exchange incorporated in the exchange−correlation potential is a major problem in UDFT approaches, the methodological conclusions of the present work are not affected by this issue. A. m-Xylylene. The vertical energy gap E*BS − Emin T is here calculated at 0.011 149 au. Let us first consider the case where both spin-decontamination factor λ and singlet-state optimized geometry Qmin S take their ideal values of 2. As δEBS and δET are here close enough (0.000 745 and 0.000 693 au, respectively), one can take a mean value for both at 0.000 719 au and use eq 6 to estimate Emin − Emin * − Emin S T from EBS T . In this way one gets min Emin − E = 0.019 422 ua (0.53 eV). Taking the explicit values S T for δEBS and δET and using eq 5 has little effect on this result min since it leads to Emin S − ET = 0.019 512 ua (0.53 eV). In this min case, the position of QS along the coordinate Q can be evaluated from eq 3, which is now

Q Smin =

δET (au)

δEBS (au)

SDFa

modelb

explicit curvec

m-xylylene

0.000 693

0.000 745

2.00 1.90 1.96

1.87 1.79 1.84

1.87 1.77 1.83

p-benzyne

0.000 869

0.000 711

2.00 1.90 1.91

2.57 2.38 2.39

2.43 1.87 1.88

Spin-decontamination factor, taken at 2, λ′, or λ from top to bottom, respectively, with average ⟨S2⟩ between QT and QBS coordinates. bSee eq 3. cSee Figures 3 and 4. Here, λ′ and λ are calculated from ⟨S2⟩ values at each point of the curves.

a

Table 2. Calculated CC Bond Lengths (Å) in m-Xylylene

2δE BS 2δE BS − δE T

CCSDa a b c d

Qmin S

leading to a value of = 1.87. Going a step further will use values of SDF different from 2, either λ′ or λ. Because ⟨S2⟩ has min not exactly the same values at Xmin T (Q=0) and XBS (Q=1), one has to take averaged values of them to evaluate λ′ and λ (for mxylylene, these changes are tiny, but we shall see later on that this is not always the case). In doing so, one gets λ′ = 1.90 and λ = 1.96. Reporting these values into eq 3 leads to values of Qmin of 1.79 and 1.84, respectively. From eq 4, the adiabatic S min energy gap Emin S − ET is now estimated at 0.018 653 au (0.51 eV) when λ′ = 1.90, and 0.019 170 au (0.52 eV) when λ = 1.96. Time has now come to check how these predictions compare to the values obtained from explicit calculations along the Q collective coordinate. To do so, we consider a simultaneous and regular variation of all internal coordinates from the triplet-state geometry Xmin T (Q=0) to that of the broken-symmetry Ms = 0 solution Xmin BS (Q=1). For each parameter of the Z-matrix, the difference between the bounds of the [0, 1] interval along Q is divided into ten regular elementary increments. Exploring the potential curve along this collective coordinate up to Q=2 and beyond therefore requires to increment more than 20 times from the starting triplet geometry. In this way we get explicit calculations for a set of regularly interpolated and extrapolated points along the collective coordinate. On each of these points, two UDFT calculations are carried out, one for the triplet state (S=1), the other one for the broken-symmetry Ms=0 solution. The corresponding results are reported in the two lower curves of Figure 2, both exhibiting clear quadratic shape. For each of these points, the spin-decontaminated singlet-state energy can be evaluated straightly from the eq 7, leading to the bold upper curve in Figure 2. Here, the spin-decontamination factor is no longer averaged but is obtained using the values of ⟨S2⟩T and ⟨S2⟩BS calculated at each point of the curve. Using 2, λ′, or λ as SDF (λ is used in Figure 2) has here little effect on the minimum of the curve, located at Q=1.87, 1.77, or 1.83 in the three cases respectively, in remarkable agreement with the results given by the above-discussed model, as summarized in Table 1. The so-extrapolated singlet-state geometrical parameters corresponding to Qmin are listed in Table 2. On the S

a

UB3LYP/6-311G**

triplet

singlet

triplet

BS Ms=0

extrapolated singlet

1.392 1.433 1.419 1.402

1.389 1.402 1.400 1.452

1.389 1.432 1.418 1.400

1.389 1.419 1.411 1.420

1.389 1.406 1.403 1.440

From EOM-SF-CCSD with UHF/6-31G* reference (ref 13b).

overall, they are in good agreement with benchmark values calculated at EOM-SF-CCSD level,13b illustrating how the whole deformation occurring from the triplet UDFT geometry to that of the BS Ms=0 solution must be continued by about the same extent. The adiabatic energy difference Emin − Emin is obtained S T directly from the minima of the triplet and spin-decontaminated singlet curves. For SDF taken at 2, λ′, and λ, one gets adiabatic gaps of 0.019 514 au (0.53 eV), 0.018 701 au (0.51 eV), and 0.019 202 au (0.052 eV), respectively, in close agreement with the values given by the above model, and exhibiting small corrections only with respect to simple brokensymmetry calculations, as illustrated in Table 4. B. p-Benzyne. As a σ diradical with a singlet ground state, p-benzyne is a good test to check whether the above-developed procedures can be likewise appropriate when the ground state is singlet. Due to this state of thing, the situation depicted in Figure 1 is now reversed. δEBS and δET still have the same definition and are here calculated at 0.000 711 and 0.000 869 au, respectively. The starting energy gap is now the vertical gap EBS * − Emin T = 0.005 064 au. From an average value of δEBS and δET, eq 6, which is now written min E Tmin − ESmin = 2(E T* − E BS ) − 4δE BS

leads to an adiabatic singlet-to-triplet gap of 0.006 986 au (0.19 eV). Using the exact values for δEBS and δET, eq 5 lead to a min modified gap of Emin = 0.006 471 au (0.18 eV). Again, T − ES 2 averaging the values of ⟨S ⟩T and ⟨S2⟩BS for Q=0 and Q=1 leads to mean values for λ′ and λ at 1.90 and 1.91, respectively. min Obtained from eq 3, Qmin S is now significantly beyond 2, (QS = min min 2.38 and 2.39, respectively), giving a ET − ES gap from eq 4 at 0.006 414 and 0.006 419 au, respectively (0.17 eV). 8230



Article

same amount as its triplet distance provides a reasonable geometry for the singlet state, to be used in a subsequent singlepoint calculation. In this way, one would obtain here a value for min Emin quite close to that given above (0.009 216 au or T − ES 0.25 eV). Because of the singlet ground state in p-benzyne, the correction brought to the simple spin-decontaminated calculations using BS Ms=0 geometry here tends to increase the adiabatic gap, whereas for m-xylylene, it tended to decrease it. As summarized in Table 4, these corrections remain weak in such small systems (less than 0.02 eV) but have a significant impact on the vertical transition energies (0.08 eV). C. Diphenalenylacetylene. Consisting of two phenalene radicals linked by an acetylene bridge, this system has been studied by Nakano et al. either at the UB3LYP level15a or by using long-range corrections to the DFT treatment (LCUB3LYP).15b Its architecture ensures an antiferromagnetic coupling according to Ovchinnikov’s rule.17 It belongs to a variety of organic architecture tailoring couplings of spin carriers to reach antiferro-, ferro-, or ferrimagnetic organizations.18 Calculations show a BS Ms=0 lowest-energy solution with spins localized on each phenalene group, the corresponding spin densities being ±0.86 electron.15b Provided that the bridge becomes cumulenic, this molecule accepts a Kékulé spin pairing into double bonds, each phenalene unit bearing in this case a bond-alternation part and a naphthalene moiety (Scheme 1, right). It is not clear therefore that the system will keep a diradicalar character. Actually, at the triplet geometry, a closed-shell solution is below the UDFT Ms=1 energy, but the BS Ms=0 solution is the lowest one. The mean value of the S2 operator is here unusually low (⟨S2⟩ = 0.79), indicating the important weight of the closed-shell component. From this BS solution at triplet-state geometry, a blind optimization leads to a geometry for which the determinant is closed-shell. The corresponding geometrical changes are rather large (up to 0.05 Å, see Table 5) and reflect bond-length variations in the direction of the cumulenic form of Scheme 1. Ignoring the spin-contamination problem and its possible geometrical impact, one could conclude that this arrangement does not present a diradicalar character. However, the spin-decontaminated energy at the triplet geometry is already below that of the closed-shell singlet state at its optimized geometry. Therefore, the closed-shell solution cannot be the optimal one. Exploring the energy curves along the collective geometrical coordinate between triplet and closed-shell singlet geometries leads to a spin-decontaminated singlet solution between these two geometries. Reoptimizing the BS Ms=0 solution at its minimum position on this coordinate actually brings little geometrical change, so that the BS Ms=0 solution can be considered to lie rather exactly on the collective coordinate linking the triplet state to the closed-shell singlet state, more precisely at 77% toward the closed-shell geometry. In other min words, the vectors Xmin → Xmin → Xmin T BS and XT closed‑shell are practically collinear. Renormalizing the collective coordinate as above with Q=0 for the triplet geometry and Q=1 for the BS Ms=0 solution leads to the energy curves plotted in Figure 5. To get the spin-decontaminated singlet curve, it is necessary, here more than ever, to use a spin-decontamination factor λ calculated at each point of the curve, because this parameter changes dramatically along the Q coordinate, as shown in the bottom curve of Figure 3 (again, if ⟨S2⟩T remains virtually constant at 2.07∼2.08 from Q=0 to Q=1, ⟨S2⟩BS varies along this interval from 0.79 to 0.42). As a consequence of this

We next proceed to the complete scanning over the entire [0−2] interval on the Q coordinate, as done in the preceding case. In contrast to it, ⟨S2⟩BS here varies substantially along the coordinate, as reflected by the λ curve reported in Figure 3, top

Figure 3. Variation of the spin-decontamination factor λ along the collective geometrical coordinate Q.

(from Q=0 to Q=2, ⟨S2⟩BS varies from 0.97 to 0.91, whereas ⟨S2⟩T remains virtually constant at 2.01). Consequently, in calculating the spin-decontaminated singlet curve from eq 7, better written here as Es = E T − λ(E T − E BS)

it is crucial to use λ values explicitly calculated at each point. In this way, one gets the bold lower curve plotted in Figure 4. Note its quadratic shape, as in BS and triplet curves, and its minimum close to Q=2, namely Qmin =1.88. Again, the soS obtained geometrical parameters are in good agreement with literature values calculated for the singlet state, as illustrated in Table 3. Using λ′ instead of λ as SDF does not change significantly this value (Qmin S =1.87), whereas using the constant factor 2 in the above equation has a dramatic impact, leading to Qmin = 2.43, a value close to that obtained above from the S model, as summarized in Table 2. From the minima of triplet and spin-decontaminated singlet curves in Figure 4, one obtains min an adiabatic energy gap Emin T − ES of 0.009 232 au (0.25 eV). Attention is drawn on the fact that, here as in the preceding case, the value of Qmin obtained from the explicitly calculated S curve occurs around Q=2, validating the model developed in the first section of the paper. Consequently, simple extrapolation of the broken-symmetry Ms=0 geometry by a 8231



Article

Figure 4. Energy curves along the collective geometrical coordinate Q for p-benzyne.

Table 3. Calculated CC Bond Lengths (Å) in p-Benzyne

ref 14c

a

UB3LYP/6-311G**

triplet

singlet

triplet

BS Ms=0

extrapolated singlet

aa b

1.370 1.398

1.358 1.419

1.377 1.407

1.370 1.421

1.363 1.435

ab b

1.383 1.406

1.374 1.424

ac b

1.387 1.403

1.383 1.411

SF-DFT (50,50)/6-31G*. (8,8)/6-31G*.

exactly coincides with the closed-shell singlet one. Because the variations of ⟨S2⟩BS and λ are smooth and regular, the junction between these two curves, not drawn in Figure 5, should behave in the same way. Listed in Table 5, the resulting interpolated bond length for the singlet state differ by about 0.01 Å from those of the BS Ms = 0 solution. Notice in this table the overall agreement between the present UDFT results and those calculated at LR-UDFT level.15b At this Qmin geometry, the diradicalar character, S reflected by ⟨S2⟩BS = 0.54, corresponds to a summed spin density of ±0.60 electron on each phenalene fragment, whereas this would reduce to ±0.50 electron for the BS Ms=0 structure at Q=1, as illustrated in Figure 6. The interpolated curve of Figure 5, differing incidentally from the extrapolations around Q=2 in the prior cases, now makes it possible to bring our corrections to the singlet−triplet energy gaps. The results are summarized in Table 4, lower part (for the obvious reasons discussed above, the line corresponding to SDF = 2 is here meaningless and must be disregarded). With respect to the broken-symmetry approach, our corrected values show an increase of the adiabatic gap by about 0.05 eV, and a decrease of the vertical gap by about 0.08 eV, the latter effect being here opposite to that observed in the benzyne diradical.

b

EOM-SF-CCSD/6-31G*.

c

CASSCF-

variation, the above-developed analytic formulas cannot be used to predict the spin-decontaminated singlet geometry Qmin S , the search of which must be numerical and proceed from interpolation of the spin-decontaminated curve between Q=0 to Q=1. From Figure 5, this leads to a value of Qmin S =0.65. As a consequence of the vanishing of ⟨S2⟩BS around Q≈1.6, both BS Ms=0 and closed-shell singlet curves merge beyond this point. There, λ = 1 so that the spin-decontaminated singlet-state curve

4. POSSIBLE GENERALIZATIONS Seeking whether this strategy might be useful for other magnetic architectures beyond interactions between two 1/2 spins, two possible extensions come to mind. The first one concerns the interaction between sites bearing spins larger than 8232



Article

Table 4. Calculated Singlet/Triplet Transition Energies (eV) broken-symmetry geometry SDF

vertical

adiabatic

vertical

adiabatic

2 λ′ λ

0.61 0.60 0.59

0.55 0.54 0.54

p-benzyne (1Ag → 3B1u)

2 λ′ λ

0.28 0.26 0.26

0.25 0.24 0.24

0.36 0.33 0.34

0.28 0.25 0.25

diphenalenylacetylene (1Ag → 3Bu)

2 λ′ λ

0.86 0.52 0.54

0.60 0.36 0.37

0.67 0.45 0.46

0.61 0.39 0.41

1

0.53 0.52 0.52

Spin decontamination factor (see the text for the definition of λ′ and λ). bDetermined for SDF = λ.

determinants may be considered, namely the solution with Ms = 1/2(m + n), which gives the energy Em+n of the upper multiplet of the AB dimer, and the solution of Ms = 1/2(m − n), where m α-electrons are centered on site A and n βelectrons are centered on site B. This solution, of energy Em−n, is not a spin eigenstate, but the corresponding energy difference enables one to fix the value of the magnetic coupling, through the relation Em+n − Em−n = 1/2 JAB mn. Knowing the value of JAB leads to the full low-energy spectrum of the AB complex. Again, optimizations of the two solutions lead to different geometries. From them one may define a collective geometry change and explore the energies of the various eigenstates of the Heisenberg Hamiltonian along this coordinate. Applying this to an antiferromagnetic interaction between two S=1 spins, and neglecting spin contamination, one may write EQ − ES = 3JAB = 3/2(E4 − E0). Again, if the curvatures of the two solutions along the Q coordinate are similar, one expects a singlet-state equilibrium geometry around QSmin = 3/2, suggesting a singlet-state relaxation from quintet geometry at 9/4 of that of the BS solution. To illustrate this procedure, we can study the larger analog of 3, ditriangulenyl acetylene, in which two triangulene units of spin S = 1 are linked by an acetylene bridge joining middle-edge carbons.19 The coupling between the two units is here antiferromagnetic, as occurring with the phenalenyl groups in 3. Similarly, quinonization of the bonds would result in a lowering (here from 1 to 1/2) of the spins carried by each unit (Scheme 2). Actually, from the Ms=2 solution to the Ms=0 one, the bond-length changes along the bridge have the right signs, but their amplitudes remain moderate (−0.04, +0.02, −0.04 Å along the bridge; altogether, only ten or so CC bond lengths vary by more than 0.01 Å). The mean value of S2 for the Ms=2 solution is slightly larger than 6, due to the spin polarization of the closed shells (⟨S2⟩ = 6.19). For the Ms=0 solution, it is significantly smaller than 2, due to the ionic components of the wave function, and the tendency to reduce the number of unpaired electrons (⟨S2⟩ = 1.67). It diminishes when going from the equilibrium geometry of the quintet (⟨S2⟩ = 1.90) to the optimal geometry of the Ms=0 solution (⟨S2⟩ = 1.67) and beyond. If these deviations are taken into account through a λ factor as proceeded above, the energy difference can be written as

Scheme 1

Table 5. Selected CC Bond Lengths (Å) in Diphenalenylacetylene

LR-UDFTa

a b c d e f a

spin-decontaminated geometryb

m-xylylene ( B2 → A1) 3

a

a

UDFT

triplet

BS Ms=0

triplet

BS Ms=0

closed-shell singlet

interpolated singlet

1.221 1.418 1.409 1.384 1.420 1.439

1.237 1.382 1.431 1.366 1.436 1.457

1.215 1.416 1.407 1.381 1.418 1.437

1.232 1.376 1.431 1.362 1.436 1.457

1.238 1.365 1.438 1.357 1.441 1.463

1.226 1.391 1.423 1.369 1.429 1.450

Reference 15b.

1/2, the magnetic interaction of which can be studied from the BS-UDFT approach. Let us consider two sites A and B, bearing m and n unpaired electrons, respectively, with maximum ground-state spin multiplicity at each center. The corresponding magnetic-coupling amplitude can be evaluated from the energies of two single determinants with m + n open shells that belong to the set of vectors on which the Heisenberg Hamiltonian HH = JABS⃗A·S⃗B is defined (i.e., the products of ground-state single determinants on each site). These energies will be identified with its diagonal terms. Two such single

EQ − ES = 3JAB = λ(EQ − E BS) 8233



Article

Figure 5. Energy curves along the collective geometrical coordinate Q for diphenalenylacetylene.

EQ − ES =

⟨S2⟩Q 2

2

⟨S ⟩Q − ⟨S ⟩BS

distributions of spins on the various sites. For instance, with M sites bearing 1/2 spins, starting from the upper multiplet of Ms=(1/2)M (all spins parallel), one will consider the M solutions where one spin has been reversed (Ms=(1/2)M − 1), and so on. From these energies it is possible to extract the (1/2)M(M + 1) intersite magnetic coupling constants, again through a mapping of these energies on the diagonal energies of the Heisenberg Hamiltonian. Typically, in a three-site problem of 1/2 spins, one will first optimize the geometries of the four solutions corresponding to the spin distributions abc, abc, ̅ ab̅c, abc.̅ The geometry changes from the Ms=3/2 solution to each of the three Ms=1/2 solutions then define three collective-motion vectors, which can be orthogonalized in the R3N−6 space if N is the total number of atoms. From these orthogonal vectors, it becomes possible to explore the energy of the various eigenstates of the Heisenberg Hamiltonian defined for each value of the coordinate along these vectors. As the number of atoms becomes large, this procedure could avoid tedious full geometry optimizations.

(EQ − E BS)

Figure 7 reports potential energy curves for the BS Ms=0 solution, the quintet state, and the so-decontaminated singlet state. Its minimum occurs at Q ≈ 1.3, a value somewhat smaller than the above crude value of 3/2. This is a consequence of the large spin contamination of the Ms=0 solution. From the singlet and quintet curves, it is straightforward to draw the triplet curve, located above the singlet, at one-third (1J) of their gap (3J). By doing so, one obtains the dashed curve in Figure 7, which has its minimum at Q ≈ 0.75, a value somewhat smaller than the ideal value of Q ≈ 1. Indeed, as the triplet state is an out-of-phase combination of two single determinants that do not interact, its equilibrium geometry is expected to be close to that of the BS Ms=0 solution. Here also, the discrepancy probably originates in the large spin contamination of the Ms=0 solution. From the decontaminated singlet minimum, the singlet-toquintet transition energies come out at 0.65 eV (vertical) and 0.48 eV (adiabatic), to be compared with the values obtained from the BS Ms=0 minimum (0.57 and 0.47 eV, respectively). The singlet-to-triplet transition energies are calculated at 0.22 eV (vertical) and 0.19 eV (adiabatic). A second possible generalization concerns systems of more than two magnetic sites. For such systems, the BS-UDFT strategy consists in finding as many solutions as independent

5. CONCLUDING REMARKS The present work has shown that spin decontamination of BSUDFT solutions may lead to significant geometry changes. If these are negligible in coordination chemistry, where magnetic orbitals are atomic d orbitals centered on transition metal ions, they may be important in conjugated organic polyradicals, where magnetic units involve both delocalized orbitals and a large number of internal coordinates. Since conjugated bridges 8234



Article

two solutions at these two geometries is sufficient to estimate the optimal geometry and the energy of the spin-decontaminated singlet. The accuracy of the prediction is easily checked by performing a calculation at the analytically predicted geometry or around it. These statements have been checked and illustrated on three different diradicalar conjugated hydrocarbons, for which the geometry changes between the singlet and the triplet states are far from being negligible. They represent three situations: (i) a triplet ground state; (ii) a singlet ground state of plain openshell character, not accepting pairing of the electrons; (iii) a system accepting such a pairing and potentially losing its diradicalar character. In the first two cases the geometry changes are larger than those predicted by simple geometry optimization of the broken-symmetry solution. They are significant and have non-negligible impact on vertical and adiabatic singlet-to-triplet energy differences. Comparisons with wave function theory-based calculations show here a rather good agreement regarding geometry changes and energetic impact. In both cases the simple analytic models predicting the geometry changes and the energy of the spin-decontaminated singlet state from unrestricted solutions have proven to be valid, as checked by recalculating the energy of the spindecontaminated state at its predicted optimal geometry. These simple models are based on a negligible variation of the spin-decontamination factor under the concerned geometry changes. This assumption is manifestly irrelevant in the third problem, where the closed-shell component of the BS solution is large and increases when the geometry changes from the BS one to the spin-decontaminated one. The above remarks should affect the debate concerning the possibly diradicalar character of long polyacenes. Calculated from DFT-optimized geometries, singlet-to-triplet energy differences can be questioned if they rely on poor determinations of the singlet geometry, either from closedshell solutions or from spin-contaminated unrestricted BS Ms = 0 single-determinantal treatments.20−22 In particular, one may notice that sophisticated wave function theory-based calculations leading to opposite conclusions rely on different DFTbased geometry optimizations.23−26 Their relevance could be examined along with the here-proposed procedure. As it is, the present strategy may be used directly for systems bridging two high-spin units (S > 1/2), such as Ni d8 ions or high-spin hydrocarbons.18,27 For systems involving more than two magnetic units as in A−B−C architectures, UDFT is frequently employed by mapping different BS-solution energies on those of Ising Hamiltonian. The principle of generalization of the present strategy to such polyradicalar systems, where several geometry changes have to be considered and combined, has been outlined and will be addressed in forthcoming work. A last comment will concern thermodynamic magnetic properties. We saw that geometry optimization results in significant variations of vertical and adiabatic differences between the singlet and triplet surfaces. As the curvatures of the potential-energy surfaces are different, this adds an extra entropic factor to the spin-degeneracy contribution. Thermodynamic magnetic properties such as magnetic susceptibility should be affected accordingly, and the use of a simple Heisenberg Hamiltonian to fit their temperature dependence should be questioned.

Figure 6. Summed spin-density in one phenalenyl fragment along the collective geometrical coordinate Q, calculated for the brokensymmetry Ms = 0 solution. The arrow indicates the spindecontaminated singlet optimal geometry.

Scheme 2

also introduce some nuclear degrees of freedom, the geometry changes from high- to low-multiplicity eigenstates a priori concern a large number of nuclear coordinates. The conception of rational (even if approximated) exploration tools of these geometry changes and of their impact on the vertical and adiabatic transition energies is certainly welcome. Although more rigorous solutions, involving derivatives of the decontamination factor with respect to nuclear coordinates are conceivable and have been proposed,9 the present work suggests extremely simple recipes to get reasonable estimates of optimized geometries of spin-decontaminated solutions and of their energies. In most cases, when the decontamination factor does not change much between the optimized geometries of the Ms=1 and Ms=0 BS solutions, calculation of the 8235



Article

Figure 7. Energy curves along the collective geometrical coordinate Q for ditriangulenylacetylene.

■

Yamaguchi, K. Chem. Phys. Lett. 2009, 483, 168. (c) Kitagawa, Y.; Saito, T.; Nakanishi, Y.; Kataoka, T.; Matsui, T.; Kawakami, T.; Okumura, M.; Yamaguchi, K. J. Phys. Chem. A 2009, 113, 15041. (d) Kitagawa, Y.; Saito, T.; Nakanishi, Y.; Kataoka, Y.; Matsui, T.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Int. J. Quantum Chem. 2010, 110, 3053. (10) Tsuchimochi, T.; Scuseria, G. E. J. Chem. Phys. 2011, 134, 064101. (11) See for instance: Kahn, O. Molecular Magnetism; VCH: Weinheim, 1993. (12) Calzado, C. J.; Cabrero, J.; Malrieu, J.-P.; Caballol, R. J. Chem. Phys. 2002, 116, 2728. (13) (a) Wenthold, P. G.; Kim, J. B.; Lineberger, W. C. J. Am. Chem. Soc. 1997, 119, 1354. (b) Wang, T.; Krylov, A. I. J. Chem. Phys. 2005, 123, 104304. (14) (a) Wenthold, P. G.; Squires, R. S.; Lineberger, W. C. J. Am. Chem. Soc. 1998, 120, 5279. (b) Sander, W. Acc. Chem. Res. 1999, 32, 669. (c) Epifanovsky, E.; Krylov, A. I. Mol. Phys. 2007, 105, 2515. (d) Wang, E. B.; Parish, C. A.; Lischka, H. J. Chem. Phys. 2008, 129, 044306. (e) For higher polyacene analogs of p-benzyne, see: Kim, H.J.; Wang, X.; Ma, J.; Cho, J.-H. Chem. Phys. Lett. 2011, 516, 141. (15) (a) Ohta, S.; Nakano, M.; Kubo, T.; Kamada, K.; Ohta, K.; Kishi, R.; Nakagawa, N.; Champagne, B.; Botek, E.; Takebe, A.; Umezaki, S.; Nate, M.; Takahashi, H.; Furukawa, S.; Morita, Y.; Nakasuji, K.; Yamaguchi, K. J. Phys. Chem. A 2007, 111, 3633. (b) Nakano, M.; Kishi, R.; Yoneda, K.; Inoue, Y.; Inui, T.; Shigeta, Y.;

AUTHOR INFORMATION

Corresponding Author

*Fax: +33 (0)5 61 55 60 65. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) For reviews, see: Neese, F. Coord. Chem. Rev. 2009, 253, 526. Bencini, A. Inorg. Chim. Acta 2008, 361, 3820. (2) Noodleman, L. J. Chem. Phys. 1981, 74, 5737. (3) Noodleman, L.; Davidson, E. R. Chem. Phys. 1986, 109, 131. (4) Noodleman, L. Inorg. Chem. 1988, 27, 1977. (5) L. Noodleman, N.; Baerends, E. J. J. Am. Chem. Soc. 1984, 106, 2316. (6) Yamaguchi, K.; Takahara, Y.; Fueno, T. In Applied Quantum Chemistry; Smith, V. H., Ed.; V. Reidel: Dordrecht, The Netherlands, 1986; p 15. (7) Soda, T.; Kitagawa, Y.; Onishi, T.; Tanako, Y.; Shigeta, Y.; Nagao, H.; Yoshika, Y.; Yamaguchi, K. Chem. Phys. Lett. 2000, 319, 223. (8) Caballol, R.; Castell, O.; Illas, F.; Moreira, I. P. R.; Malrieu, J. P. J. Phys. Chem. A 1997, 101, 7860. (9) (a) Kitagawa, Y.; Saito, T.; Ito, M.; Shoji, M.; Koizimi, K.; Yamanaka, S.; Kawakami, T.; Okomura, M.; Yamaguchi, K. Chem. Phys. Lett. 2007, 442, 445. (b) Saito, T.; Nishihara, S.; Kataoka, Y.; Nakanishi, Y.; Matsui, T.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; 8236



Article

Kubo, T.; Champagne, B. J. Phys. Chem. A 2011, 115, 8767. (c) For differently-coupled two-phenalene systems, see: Nakano, M.; Nakagawa, N.; Kishi, R.; Ohta, S.; Nate, M.; Takahashi, H.; Kubo, T.; Kamada, K.; Ohta, K.; Champagne, B.; Botek, E.; Morita, Y.; Nakasuji, K.; Yamaguchi, K. J. Phys. Chem. A 2007, 111, 9102. (d) For staggered dimeric pairs of phenalenyl radicals, see: Takano, Y.; Taniguchi, T.; Isobe, H.; Kubo, T.; Morita, Y.; Yamamoto, K.; Nakasuji, K.; Takui, T.; Yamaguchi, K. J. Am. Chem. Soc. 2002, 124, 11122. (16) Frisch, M. J.; et al. Gaussian 03, revision B.05; Gaussian, Inc.: Wallingford, CT, 2004. (17) Ovchinnikov, A. A. Theor. Chim. Acta 1978, 47, 297. Also known in solid-state physics as the Lieb’s theorem: Lieb, E. H. Phys. Rev. Lett. 1989, 62, 1201. (18) Trinquier, G.; Suaud, N.; Guihéry, N.; Malrieu, J.-P. Chem. Phys. Chem. 2011, 12, 3020. (19) For the weak energy-impact nonplanar distortion of this system, see ref 18, p 3029. (20) Bendikov, M.; Duong, H. M.; Starkey, K.; Houk, K. N.; Carter, E. A.; Wudl, F. J. J. Am. Chem. Soc. 2004, 126, 7416. (21) Jiang, D.; Dai, S. J. Phys. Chem. A 2008, 112, 332. (22) Dos Santos, M. C. Phys. Rev. B 2006, 74, 045426. (23) Hachmann, J.; Dorando, J. J.; Avilès, M.; Chan, G. K.-L. J. Chem. Phys. 2007, 127, 134309. (24) Hajgato, B.; Szieberth, D.; Geerlings, P.; De Proft, F.; Deleuze, M. S. J. Chem. Phys. 2009, 131, 224321. (25) Hajgato, B.; Huzak, M.; Deleuze, M. S. J. Phys. Chem. A 2011, 115, 9282. (26) Huzak, M.; Deleuze, M. S.; Hajgato, B. J. Chem. Phys. 2011, 135, 104704. (27) Trinquier, G.; Suaud, N; Malrieu, J.-P. Chem. Eur. J. 2010, 16, 8762.

8237


A Recipe for Geometry Optimization of Diradicalar Singlet States from

Recommend Documents