Characterizing Bonding Patterns in Diradicals and ... - ACS Publications

Subscriber access provided by READING UNIV

Article

Characterizing bonding patterns in diradicals and triradicals by density-based wave function analysis: A uniform approach Natalie Orms, Dirk Rehn, Andreas Dreuw, and Anna I. Krylov J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01012 • Publication Date (Web): 21 Dec 2017 Downloaded from http://pubs.acs.org on December 26, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Characterizing Bonding Patterns in Diradicals and Triradicals by Density-Based Wave Function Analysis: A Uniform Approach Natalie Ormsa , Dirk Rehnb , Andreas Dreuwb , Anna I. Krylova† a Department

of Chemistry, University of Southern California, Los Angeles, California 90089-0482, U.S.A. b

Ruprecht-Karls University of Heidelberg and

Interdisciplinary Center for Scientific Computing, Im Neuenheimer Feld 205A, D-69120, Heidelberg, Germany

December 15, 2017 Abstract Density-based wave function analysis enables unambiguous comparisons of electronic structure computed by different methods and removes ambiguity of orbital choices. We use this tool to investigate the performance of different spin-flip methods for several prototypical diradicals and triradicals. In contrast to previous calibration studies that focused on energy gaps between high and low spin-states, we focus on the properties of the underlying wave functions, such as the number of effectively unpaired electrons. Comparison of different density functional and wave function theory results provides insight into the performance of the different methods when applied to strongly correlated systems such as polyradicals. We show that canonical molecular orbitals for species like large coppercontaining diradicals fail to correctly represent the underlying electronic structure due to highly non-Koopmans character, while density-based analysis of the same wave function delivers a clear picture of bonding pattern.

1

Introduction

Chemists define diradicals and triradicals as species with two and three unpaired electrons.1–4 This bonding pattern arises when two (or three) electrons are distributed in two (three) nearly degenerate molecular orbitals. Fig. 1 shows all resulting configurations with positive spinprojections. For same-spin electrons, there are only two possible arrangements: Ms ± 1 (configuration A(i) in Fig. 1 has Ms =1) and Ms ± 3/2 (configuration B(i) in Fig. 1 has Ms =3/2). For the states with lower spin-projections, more configurations can be generated, as illustrated in The energy gaps, relative state ordering, and relative weights of the Slater determinants (i.e., coefficients λ in Fig. 1) depend on the nature and energy separation of the frontier molecular orbitals (MOs). The character of the MOs also determines the character of the wave function,

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(i)

+

-

(ii)

(iii)

!

Page 2 of 25

+

-!

(iv)

(v)

(A) + (i)

!’

+

2

-

(iii)

(ii)

+

- !”

(vi)

(vii)

!”

-

- !’

(iv)

(v)

+

- !’’’

(viii)

(ix)

!’’’

+ (x)

(B) Figure 1: Wave functions of di- (A) and triradicals (B) that are eigenfunctions of S 2 . In both panels, wave function (i) corresponds to the reference-state wave function: the high-spin Ms = 1 triplet state of a diradical or the Ms = 3/2 quartet state of a triradical. Wave functions (ii)-(v) in A and (ii)-(x) in B correspond to the low-spin states: Ms = 0 singlets and triplets, or Ms = 1/2 doublets and quartets. For simplicity, only configurations with positive spin-projection are shown.

e.g., whether it has predominantly covalent (i.e., two electrons residing on different parts of a molecule) or ionic character. To illustrate this point, consider a simple diradical, such as a molecule with a broken bond (i.e., H2 at dissociation limit or twisted ethylene). In this case, the two MOs are a bonding and antibonding pair (such as σ and σ ∗ orbitals in stretched H2 or π and π ∗ in twisted ethylene). When the two MOs are exactly degenerate, λ=1 in wave functions (iv) and (v) in Fig. 1A. Wave function (iii) corresponds to the covalent diradical singlet state (H· H· ), which is exactly degenerate with the triplet state, whereas (v) represents a purely ionic (charge-resonance) state (i.e., H− H+ ± H+ H− ). This situation is often described as a perfect diradical. If the gap between bonding and antibonding orbitals is large, λ becomes small and the ground state can be represented by a single Slater determinant; the wave function in this case is a mixture of covalent and charge-resonance configurations and the molecule can be described as a closed-shell species. In between there is a continuum of intermediate bonding patterns, which are sometimes described as diradicaloids.2 Detailed analysis of different types of diradical electronic structure can be found in classic papers;1, 2 for a review of triradicals, one can consult Ref. 3. The wave functions composed by Slater determinants in which two electrons occupy the same MO, e.g., (iv) and (v) in (A) and (v)-(x) in (B) are commonly referred to as of “closedshell” type (although they can correspond to completely unpaired electrons, as in the case of 2


Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A(v) with λ=1), in contrast to “open-shell” wave functions (i)-(iii) in (A) and (i)-(iv) in (B) (which can correspond to a purely ionic and, therefore, closed-shell pattern). Here we follow this accepted terminology for consistency with other studies.4 An important difference between the two types of open-shell wave functions is that in A(v) the coefficients λ depend on the orbital degeneracy whereas in A(ii)-(iii) the relative weights of the contributing Slater determinants are determined by spin symmetry alone.1 In model pedagogical examples of two-electrons-in-two orbitals or three-electrons-in-threeorbitals, the number of unpaired electrons can be unambiguously determined from the wave functions. For example, in diradicals, (i) and (ii) correspond to the two unpaired electrons, whereas the relative weight of purely covalent configuration in (v) is determined by λ, i.e., λ=1 corresponds to the 2 unpaired electrons. However, this simple picture does not apply for realistic many-body wave functions. Even in the simplest case of two electrons, the dynamic correlation and the arbitrariness in orbital choice result in multi-configurational wave functions whose character cannot be easily assigned on the basis of just 2 leading configurations. In many-electron systems the contributions of electrons occupying lower MOs give rise to throughbond interactions between the radical centers, further complicating the wave function analysis. Although the number of effectively unpaired electrons is not an observable physical property, it is related to the bonding pattern, which, in turn, manifests itself in concrete structural, spectroscopic, and thermochemical quantities.3–5 Thus, the ability to assign an effective number of unpaired electrons to a particular electronic wave function is valuable for chemical insight, akin to other methods of wave function analysis.6, 7 Several solutions towards this goal have been proposed.8–10 As with many other wave function analysis tools,6, 7, 11–14 they exploit the concept of natural orbitals.15 Natural orbitals (ψi ) are eigenstates of the one-particle reduced density matrix: X ψi (x) = cji φj (x), (1) j

γpq

γc = cn, ≡ hΨ|p† q|Ψi

(2) (3)

where pˆ† and qˆ are the creation and annihilation operators corresponding to φp and φq MOs. The respective eigenvalues (ni ) are non-negative, add up to the total number of electrons, and can be interpreted as orbital occupations. Natural orbitals afford the most compact representation of the electron density X X ρ(x) = γpq φp (x)φq (x) = ni ψi2 (x) (4) pq

i

and reflect multi-configurational character of wave functions. For example, for a single Slater determinant, natural occupations are two (for occupied orbitals) and zero (for virtuals). For a twodeterminantal wave function (v) with λ=1, natural occupations are 1. For multi-configurational wave functions, the partial occupations of natural orbitals can be used to derive an effective number of the unpaired electrons. For example, for the perfectly diradical wave functions, such as A(i)-(iii) and (iv) with λ=1 in Fig. 1, n1 = n2 =1 and the number of unpaired electrons is 2, whereas for the closed-shell case, A(iv) with λ=0, n1 =2, n2 =0 and the number of unpaired electrons is 0. Several ways to compute an effective number of unpaired electrons from the one-particle density matrices have been proposed.8–10 In this work, we make use of the two indices, nu and 3



nu,nl , proposed by Head-Gordon9 as an extension of work by Yamaguchi et al.:8 X nu = min(¯ ni , 2 − n ¯i)

Page 4 of 25

(5)

i

nu,nl =

X

¯ i )2 n ¯ 2i (2 − n

(6)

i

In both equations, the sum runs over all natural orbitals and the contributions of the doubly occupied and unoccupied orbitals are exactly zero. As one can easily verify, both formulas yield correct answers for the model examples from Fig. 1, i.e., 2 for A(i)-(iii) and (iv) with λ=1, 0 for A(iv) with λ=0, and so on. The two expressions differ by how they account for partially occupied orbitals. Numerical experimentation has shown that for many-electron wave functions, Eq. (6) consistently gives physically meaningful values, whereas nu often produces the number of unpaired electrons which is too high (as compared to chemical intuition). This happens because the nu index does not suppress dynamic correlation contributions (which come from a large number of small natural occupation numbers) to the total number of unpaired electrons, whereas nu,nl emphasizes radical character at ni values near one and closed-shell character for ni values close to zero and two. Below we refer to nu,nl as Head-Gordon’s index. For the two-electrons-in-two-orbitals case, the coefficient λ and nu,nl are related by the following expression:5 32λ4 . (7) nu,nl = (1 + λ2 )4 As illustrated by the numerical examples in Ref. 5, nu,nl provides a much more robust and reliable measure of the diradical character than λ. Furthermore, unlike λ’s, nu,nl does not depend on orbitals used in the correlated calculation (unrestricted or restricted open-shell Hartree-Fock, Kohn-Sham orbitals, etc). In addition to enabling calculations of the number of effectively unpaired electrons, natural orbitals afford visualization of the frontier orbitals that are associated with correlated manyelectron wave functions, thus departing from canonical Hartree-Fock MOs. Comparing natural occupations of the frontier orbitals with the computed number of effectively unpaired electrons informs us of how well the frontier orbital picture (e.g., two-electrons-in-two-orbitals representation of diradicalas) represents the reality. In the idealized diradical case, only two frontier natural orbitals contribute to Eq. (6). Because the one-particle density matrix is a reduced quantity, which can be computed for any wave function (or for Kohn-Sham DFT and TD-DFT states), the analyses based on density matrices afford comparisons of bonding patterns computed by different methods and are also orbital-invariant. The utility of state or transition density matrices as well as difference density matrices are well recognized in quantum chemistry.6, 7, 11–14, 16–21 Here we apply state-density based analysis tools to characterize the electronic structure of several prototypical diradicals and triradicals. Polyradicals play important roles in fields as diverse as photochemistry,22 atmospheric chemistry,23 and molecular magnetism.24 Depending on the type of frontier MOs, they can be described as all-σ or all-π, σ − π, or spatially separated atom-centric. The degree of interaction and (nominal) bonding between unpaired electrons differs for each type and each molecular species. We consider the following diradicals: methylene (CH2 , same-center diradical); ortho-, meta-, and para-benzyne (σσ); 1-(2-dehydroisopropyl)-4-dehydrobenzyne (σπ), wherein radical 4


Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


electrons are localized on the tertiary carbon of the propane substituent and the para position of the benzene ring; binuclear copper complexes CUAQAC02 and CITLAT (spatially separated atom-centric). Model triradicals include 1,3,5- and 1,2,4-tridehydrobenzyne (all-σ), and 5-dehydro- and 2-dehydro-meta-xylylene (σ −π). Together, these molecules constitute a diverse benchmark set, representative of various bonding patterns one encounters in open-shell species. In this paper we employ several variants of spin-flip (SF) methods.25–29 The SF approach affords robust and accurate description of diradicals and triradicals within a black-box singlereference formalism. The performance of different variants of the SF method has been extensively benchmarked, focusing on the energy gaps between the electronic states and sometimes their structures.27–37 The most insightful comparisons are based on photoelectron spectra, which provide information about the electronic energies and the structures.34, 35 These studies illustrated that the SF-CCSD method provides reliable energy gaps, with errors close to 1 kcal/mol, as well as reliable photoelectron spectra.34, 35 Subchemical accuracy can be achieved by including the effect of higher excitations perturbatively.38 Within SF-DFT, the best performance is delivered by B5050LYP (in collinear formulation,27 where a large fraction of exact exchange is needed), and by PBE50 (within non-collinear formulation).28 However, no analysis of the underlying wave functions has been performed, with an exception of several recent studies.35, 39 One interesting question is how the level of correlation (or the functional, in SF-DFT methods) affects the effective number of unpaired electrons. This question is addressed here, by means of the analysis of the density matrices of states obtained by the different post-HartreeFock and TDDFT methods. The structure of the paper is as follows. In the next section, we describe theoretical methods and outline computational protocols. We then present our results, beginning with a comparison of the equilibrium geometries computed by Kohn-Sham and SF-TDDFT. We use meta-benzyne to illustrate the failure of standard Kohn-Sham DFT to describe singlet-state structures of strongly correlated systems. We proceed to discuss wave function properties—in particular, nu and nu,nl —using the example of H2 along the dissociation coordinate. In Section 3.3, we analyze the character of EOM-SF-CCSD wave functions of methylene, the di- and tridehydrobenzynes, and the triradical xylylenes. We then compare these results with the analysis of SF-TDDFT and SF-ADC(2) wave functions. Included in the SF-TDDFT and SF-ADC(2) results section are the CUAQAC02 and CITLAT bi-copper complexes. We conclude by summarizing relative performance of the different approaches.

2

Theoretical methods and computational details

As illustrated in Fig. 1, diradical character results in multi-configurational wave functions. This multi-configurational character arises due to electronic near-degeneracies of the frontier molecular orbitals. Standard electronic structure methods,40 which are based on the hierarchical improvements of a single-determinantal Hartree-Fock wave function, fail in situations where more than one Slater determinant has a large contribution. The Kohn-Sham DFT also breaks down in this case. Such strongly correlated systems are sometimes treated by multi-reference methods based on a multi-configurational SCF ansatz and separate description of static and dynamic correlation. Here we employ an alternative approach, the SF method,25, 41 which allows one to describe multi-configurational wave functions of the diradical and triradical types in a simple single-reference framework. 5



Page 6 of 25

The SF approach is based on the observation that high-spin states, such as Ms =1 triplet or Ms =3/2 quartet, can be well described by a single determinant and that the low-spin states (singlets and doublets) are formally single excitations from the respective high-spin reference states. Thus, in the SF approach, a high-spin state is used as the reference state and the target manifold of states is generated by applying a linear spin-flipping excitation operator to the reference state: t ˆ Ψs,t Ms =0 = RMs =−1 ΨMs =1 , ˆ Ms =−1 Ψq Ψd,q =R . Ms =1/2

Ms =3/2

(8) (9)

For diradicals, a triplet reference (Ms =1 configuration (i) in panel A of Fig. 1) is used, and a quartet reference (Ms =3/2 configuration (i) in panel B of Fig. 1) is used for triradicals. Using different approaches to describe correlation in the reference state gives rise to different SF methods. For example, applying this strategy to a Kohn-Sham determinant leads to SF-TDDFT27, 28 ˆ includes only single excitations that flip the spin of an electron). (in this case, the operator R In wave-function methods, one can use an uncorrelated reference (Hartree-Fock), which gives rise to SF-CIS,25, 42 or a correlated one, such as MP2 or CCSD. The accuracy of SF calculations can be systematically improved (up to the exact result) by increasing the level of correlation. In this paper, we consider two wave function approaches: one based on coupled-cluster theory (EOM-SF-CCSD or SF-CCSD)26, 30 and one based on the algebraic-diagrammatic construction (ADC) scheme for the polarization propagator43 (SF-ADC).29, 36 In EOM-SF-CCSD, the operaˆ includes single and double excitations that flip the spin of an electron. In the SF-ADC(2) tor R method,29 both single and double excitations are included, but the doubly excited determinants are treated only in zeroth order of perturbation theory. The mathematical structure of the 2nd order correction is similar to the one of the SF-CIS(D∞ ) method,44 however, in CIS(D∞ ) it is evaluated perturbatively only once using CIS wave functions, whereas in ADC(2) it is included already within the iterative solution for the ADC vectors, similarly to the CC2 method.45 For selected examples, we also present results for SF-TDDFT.27, 28 As one can see, all configurations that appear in electronic states shown in Figure 1 are treated on an equal footing within the SF formalism. In diradicals, the target-state manifold comprises the singlets and the Ms =0 component of the triplet state. Likewise, in triradicals the target-state manifold comprises the doublets and the low-spin (Ms =1/2) component of the quartet state. Because of the balanced description of all target states, the energy gaps between the target SF states are more accurate than energy gaps between the reference and the target states. SF methods can be employed to optimize the geometries of the target states using analytic gradients.26–29 SF wave functions can be used to construct one-particle density matrices that can then be analyzed using natural orbitals and wave function analysis tools.

2.1

Computational Details

We performed geometry optimizations of the high-spin triplet and quartet states by DFT/B5050LYP/ccpVTZ. To obtain the structures of the lowest closed-shell singlet or doublet states, we used SF-TDDFT/B5050LYP/cc-pVTZ. For methylene, we used the FCI/TZ2P structures.46 These structures are shown in Fig. 2. In calculations of the CUAQAC02 and CITLAT binuclear copper diradicals, we used X-ray crystal structures from the Cambridge Crystal Structure Database.47, 48 The perchlorate counterion from the crystal structure of CITLAT was omitted. These structures are shown in SI. 6


Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To assess the effect of the method on the structure, for selected systems, we performed additional geometry optimizations of the singlet and triplet (diradicals) or doublet and quartet (triradicals) Kohn-Sham references. Relevant Cartesian coordinates are provided in the supplemental information (SI). For all benchmark systems, we present vertical energy gaps, computed as energy differences of the target SF states at the ground-state equilibrium geometry (optimized by SF-TDDFT). Energy gaps are defined as follows: ∆E = Elow−spin − Ehigh−spin ,

(10)

i.e., negative ∆E corresponds to the low-spin ground states (singlets or doublets); note that Ehigh−spin is the energy of the Ms =0 triplet or Ms =1/2 quartet state and not the energy of the reference state. Natural orbitals and Head-Gordon’s indices of different states are also computed at the respective equilibrium geometries. We used the cc-pVTZ basis in all calculations. For selected examples, we compare the following levels of theory: EOM-SF-CCSD, SFADC(2)-s, and SF-TDDFT. In SF-TDDFT we employed the following functionals: PBE50 (50% PBE49 and 50% Hartree-Fock exchange with 100% PBE correlation), B5050LYP (50% Hartree-Fock + 8% Slater + 42% Becke50 for exchange, with 19% VWN + 81% LYP51 for correlation),27 B97,52 and LDA (Slater exchange with VWN correlation) functionals. We used the collinear formulation with B5050LYP27 and non-collinear formulation28, 53, 54 with all other functionals. We used the libwfa module12, 13 of Q-Chem to compute and visualize natural orbitals and Head-Gordon’s indices. We report nu,nul , Eq. (6), and natural occupations of the frontier orbitals. Because the present implementation of the SF methods are not spin-adapted, the SF states show some (usually small) spin-contamination even if ROHF references are employed. Consequently, α and β frontier orbitals (and the respective occupations) are slightly different. Below we report spin-average natural occupation, n ¯i: n ¯ = |nα + nβ | .

(11)

We also report the difference between the α and β natural occupations (∆ = |nα − nβ |), which provides an additional measure of spin-contamination. In all figures, we show only orbitals of the triplet states of diradicals, because the shapes of frontier natural orbitals for singlet and triplet states are indistinguishable for all systems considered in this study. We only show αorbitals, as the shapes of paired α and β natural orbitals are the same. For triradicals, we show α orbitals for the doublet and quartet states. The Q-Chem electronic structure package was used for all calculations,55, 56 and orbitals were rendered using Jmol.

3 3.1

Results and discussion Equilibrium Geometries

The extent of open-shell character has concrete structural implications,3, 4, 57 e.g., in molecules with large diradical character the structures of triplet and singlet states are rather similar, in contrast to closed-shell systems. Consequently, the choice of electronic structure method is important for obtaining accurate structures. Figure 2 shows the optimized geometries for the 7



Page 8 of 25

1.369 1.360 1.363 1.371

1.390 1.391

1.075 1.070

1.382 1.238 118.4° 110.8°

121.1° 126.9°

1.377 1.394

120.5° 122.3°

1.075 1.072

1.074 1.075

1.366 1.359

1.073 1.073 121.3° 113.8°

1.438 1.357

1.075 1.077

1.365 1.351 125.0° 139.6°

114.8° 94.9° 116.9° 116.0°

A 1.368 1.379

1.368 1.348

1.392 1.416

1.372 1.349

1.073 1.069

F

1.075 1.074

1.074 1.074

1.413 1.404

1.384 1.245 120.3 ° 129.1 ° 119.2 ° 110.8 ° 116.6 ° 118.7 °

1.362 1.377

117.2° 117.4° 1.411 1.428

133.3°

121.4° 121.2°

D

C

1.365 1.374 1.076 1.073 1.396 1.396

1.073 1.073

118.7° 1.492 119.1° 1.488

121.9 ° 126.4 °

121.2° 119.7° 1.075 1.071

118.1° 112.1°

126.8° 125.1° 116.6° 117.5°

1.074 1.074

121.6° 117.4° 121.6° 116.9°

1.388 1.386

B

123.7° 140.0° 115.2° 96.0°

1.377 1.383

124.7° 125.4°

125.9 ° 129.7 °

127.9 ° 1.072 128.4 ° 1.071

1.390 1.383 1.385 1.391

1.367 1.376 1.077 1.071

120.9 ° 121.0 °

1.372 116.2 ° 1.371 105.2 °

114.3 ° 114.0 °

1.427 1.425 1.379 1.376

1.075 1.075 1.075 1.074

121.8 ° 121.6 °

H

G

E 1.075 1.075

123.0 ° 123.1 °

1.406 1.408 122.2° 122.2°

1.388 1.385

117.7 ° 117.6 °

117.2 ° 116.9 ° 1.074 1.074 127.2 ° 127.8 °

1.428 1.433

1.357 1.354

I

Figure 2: Structures of the four benzyne diradicals (A-D), methylene (E), and four triradicals (F-I) included in this study. Numbers in italics correspond to the optimized lowest highspin state (Ms = 1 triplets for A through E and Ms = 3/2 quartets for F through I), while underlined numbers correspond to the optimized lowest energy singlet (diradicals) or doublet (triradicals) state computed by SF-TDDFT. Additional structural data for compounds D and H is presented in the SI. Compounds will be referred to by letter designation in the remainder of this work.

two lowest spin-states of each of the di- and triradicals included in this study (see Section 2.1 for details). For D and H, we also show the structures computed by the standard Kohn-Sham DFT for the low-spin states (Fig. S2 in SI). For species with well-separated radical centers, like p-benzyne and 1-isopropyl-4-benzyne, the SF optimized geometry of the singlet or doublet state closely resembles the optimized geometry of the respective high-spin state, indicative of weakly interacting unpaired electrons. For species with modest di- or triradical character, the singlet or doublet structures are markedly different from the high-spin geometries, due to the bonding interactions between the unpaired electrons.57 The structures of triplet or quartet states can be accurately computed by using standard Kohn-Sham DFT for high-spin states (these structures are very close to the structures computed for the corresponding SF-DFT states). In cases of very weak diradical/triradical character, the structures of the low-spin states can also be computed by regular Kohn-Sham DFT, however, in cases of relatively strong di-/triradical character, one needs to employ the SF-DFT approach, in order to correctly capture multi-configurational character of the underlying wave functions. The most notable example is meta-benzyne shown in Fig. 3, where the angle of separation 8


1.073 1.073

120.9° 120.8°

Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


between the radical carbons in the SF optimized singlet geometry is roughly 20 degrees less than in the structure optimized by regular closed-shell Kohn-Sham DFT. 69.3°

94.9°

1.518 Å

Kohn-Sham DFT Singlet nu,nl = 0.03

114.8°

1.987 Å

2.305 Å

SF-TDDFT Singlet

Kohn-Sham DFT Triplet

nu,nl = 0.26

nu,nl = 0.90

Figure 3: Optimized geometries of the lowest singlet state computed by regular KohnSham DFT (left) and SF-TDDFT (center), and high-spin triplet state of meta-benzyne (B5050LYP/cc-pVTZ). Distance and angle of separation between radical sites are shown. HeadGordon’s indices for the 3 structures are computed using SF-CCSD wave functions for the lowest singlet state. 1

Using a method that fails to correctly describe open-shell character can lead to catastrophically wrong structures.4 One well studied example is meta-benzyne.58 In this system, single-reference methods such as regular Kohn-Sham DFT and CCSD underestimate the diradical character and yield (incorrect) bicyclic structures, whereas methods that do not include dynamic correlation (valence bond, CASSCF) exaggerate the distance between the radical centers. The structure used in wave function analysis has a strong effect on the bonding pattern. That is, using a wrong structure will produce an incorrect number of unpaired electrons, even if an accurate electronic structure method is used for the wave function analysis. We illustrate this point using meta-benzene. Figure 3 shows optimized structures for meta-benzyne. The rightmost structure is computed by B5050LYP for the high-spin triplet state in which the two electrons are completely unpaired. Because the distance between the two radical centers is large, the singlet-state wave function shows considerable diradical character as indicated by the relatively large value of nu,nl (0.9). The central structure is computed by SF-DFT/B5050LYP for the lowest singlet state. As one can see, the distance between the two radical centers is shorter by about 0.3 ˚ A, as compared to the triplet-state structure, due to a partial bond formed by the two electrons. At this geometry, the singlet-state wave function has moderate diradical character (nu,nl =0.3). The leftmost structure is computed using regular restricted Kohn-Sham DFT/B5050LYP. Because this approach is not capable of describing diradical character, the optimized structure is bicyclic, with a short distance between the two radical centers (this structure is similar to the CCSD structure reported by Crawford and co-workers58 ). The singlet-state SF-CCSD wave function computed at this geometry shows nearly perfect closed-shell character (nu,nl =0.03). This example illustrates the importance of performing calculations at the nuclear geometry that corresponds to the correct electronic configuration of a molecule (i.e., in the case of meta-benzyne, of a singlet with moderate diradical character).

9



# of unpaired electrons

2 1.5 1 0.5

nu (CCSD) CCSDnu

nu,nl (CCSD) CCSDnu,nl

|nu – nu,nl| (CCSD) CCSDnu nu,nl

nu (B5050LYP) DFTnu

nu,nl (B5050LYP) DFTnunl

|nu – nu,nl| DFTnu-nunl (B5050LYP)

0 0.5

1.5

2.5

3.5

4.5

5.2

0.8

5.1

0.7 0.6

!E (eV)

5.0

0.5

4.9

0.4

∆E ∆E

4.8

Series1

4.7

0.3

∆nu,nl

∆nu,nl

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 25

0.2

Series2

4.6

0.1

4.5

0 0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

R (Å) P Figure 4: nu , nu,nl , and |nu - nu,nl | (top), and EB5050LY P - ECCSD (∆E) and |nB5050LY - nCCSD u,nl u,nl | (∆n, bottom) computed for the lowest singlet state of H2 along the bond-stretching coordinate using EOM-SF-CCSD and SF-TDDFT/B5050LYP with the cc-pVTZ basis set.

3.2

Head-Gordon’s Indices Along the H2 Dissociation Curve

Before proceeding to wave-function analysis in diradicals and triradicals, let us consider a model example, for which CCSD is exact, dissociation curve of the dihydrogen molecule. Figure 4 shows nu , nu,nl , |nu − nu,nl |, and the difference in total B5050LYP and CCSD state energies as a function of internuclear distances R for the lowest singlet state of H2 . Tabulated raw data is provided in the SI. As expected, around equilibrium the number of unpaired electrons is small. As the internuclear distance R increases, diradical character increases, reaching 2 at the dissociation limit. At the equilibrium and at the dissociation limit, SF-TDDFT agrees well with SF-CCSD (exact result). However, at intermediate distances, SF-TDDFT underestimates the number of unpaired electrons. Compare, for example, the two blue curves, which show the respective nu,nl s. At 2 ˚ A, the exact wave function has 1 unpaired electron, whereas the SF-TDDFT wave function has only 0.25. Only 0.5 ˚ A further SF-DFT develops open-shell character yielding nu,nl =1. This lag is observed for nu and nu,nl . Strict SF-ADC(2) does not give exact excited states due to the zeroth-order treatment of the doubly excited determinants, but improving the description of the doubly excited configurations to first order, as in SF-ADC(2)-x and SF-ADC(3), yields exact states for H2 , as SF-CCSD (note that only the energy differences, such as singlet-triplet gaps are exact, but not the total state energies). Regardless of the level of theory, at the dissociation limit the difference between nu and nu,nl is small. However, at small and intermediate distances we observe noticeable discrepancy 10


Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


between nu and nu,nl . Compared to nu,nl , around the equilibrium nu overestimates diradical character, and at longer distances it underestimates it. As expected, maximum deviation occurs around the natural orbital occupation numbers of 0.25 and 1.75 electrons, where the quadratic nature of the expression for nu,nl suppresses and enhances radical character, respectively.9 The point along the dissociation curve where the nu and nu,nl curves intersect at nu = nu,nl ≈1 depends on the level of theory. This large discrepancy between the two quantities is only observed for the singlet state — nu and nu,nl equal exactly 2 for the high-spin reference and low-spin triplet states at every point along the dissociation curve from 0.74-5.00 ˚ A (raw data is provided in SI). The lower panel of Figure 4 shows the difference between the SF-CCSD and SF-DFT total energies and the respective nu,nl . At large internuclear separation (> 3˚ A), the two methods yield identical nu,nl and the respective potential energy curves are nearly parallel. However, at shorter distances (less than 3 ˚ A), where we observe large discrepancies between the number of unpaired electrons computed by SF-CCSD and SF-DFT, we also observe large non-parallelity errors in the SF-DFT potential energy curve. This example illustrates that the errors in SFDFT energies originate from errors in underlying densities.

3.3

EOM-SF-CCSD energies and wave function character in diradicals and triradicals

Table 1 shows energy separations, nu,nl , and hS 2 i for the lowest states of di- and triradical benzynes and methylene. nu and |nα − nβ | for each state are provided in the SI. All states considered suffer very little from spin-contamination. For high-spin states, nu,nl is very close to the ideal values of 2 unpaired electrons for triplet states and 3 unpaired electrons for quartet states. In contrast, nu,nl for singlets and doublets depends very much on the nature of the di- or triradical ground state. For species with a singlet ground state, nu,nl of the ground state ranges from closed-shell values close to zero (for singlets with modest radical character) to values close to two (for open-shell singlets and strong diradical character). The same is true for triradicals with doublet ground states; namely, closed-shell doublets like those observed in F and G have nu,nl values close to 1, while open-shell doublets such as the ground state of I have nu,nl values close to 3. We observe that radical character (nu,nl ) increases and ∆E becomes more positive as the distance between the radical carbons increases from A through D. Natural orbitals of the lowest singlet/doublet and triplet/quartet states of methylene and the di- and triradical benzynes are shown in Figure 5. With the exception of methylene, D, H, and I, all di- and triradical frontier orbitals are of the σ type, consistent with the molecular orbital patterns reported by Krylov and Cristian.57 The natural orbitals in D, H, and I are of σπ type. The values of n ¯ show spin-averaged occupancy of each spatial orbital. The trends in orbital occupancy are consistent with nu,nl : increased radical character is ascribed to ground states with natural orbitals n ¯ -values close to 1. Values in parentheses (∆n) indicate spin imbalance in orbital occupancy arising due to spin-incompleteness of the underlying wave function. For diradicals, ∆n-values are close to the ideal value of zero. Doublet tridehydrobenzynes also exhibit ideal ∆n-values, with only one predominantly singly occupied natural orbital hosting the odd electron that gives rise to the doublet. Quartet states of the tridehydrobenzynes, and open-shell doublet and quartet states of the xylylene triradicals, feature natural orbitals that are predominantly singly occupied (¯ n close to 1), but suffer from some spin-imbalance. This 11



imbalance appears to be minor and does not strongly impact other properties of these states, like the hS 2 i values or the sign and the magnitude of ∆E. For all systems, we observe that the difference between nu,nl and nu is much smaller for the low-spin triplet states than for the singlet states. This is partially because in the limiting case of 2 unpaired electrons the difference between the two indices is expected to vanish. For example, the two indices give identical results for the singlet state of dihydrogen at the dissociation limit and identical results for the triplet states at all distances. Yet, for many-electron electron examples, even for triplet states nu assumes values that are slightly larger than 2. This can be attributed to the fact that nu,nl suppresses contributions from the dynamic correlation. Thus, the second reason for much smaller discrepancy between the two indices for the triplet states relative to the singlets can be attributed to the smaller dynamic correlation, which is characteristic of the triplet wave functions.25 Interestingly, the high-spin reference states show larger difference between nu,nl and nu than the respective low-spin components of the same multiplicity. Table 1: EOM-SF-CCSD energy gaps (eV) and wave function properties of lowest two states (singlets and Ms = 0 triplets for diradicals, Ms = 1/2 doublets and quartets for triradicals). Molecule

∆E

CH2 A B C D F G H I

0.94 -2.46 -1.90 -0.24 0.24 -2.21 -2.00 0.41 -0.14

3.4

nu,nl Singlet/Doublet Triplet/Quartet 0.25 2.00 0.16 2.00 0.26 2.00 1.45 2.00 2.00 2.01 1.28 3.06 1.16 3.02 3.01 3.02 3.01 3.01

hS 2 i Singlet/Doublet Triplet/Quartet 0.00 2.00 0.00 2.00 0.01 2.01 0.01 2.01 0.01 2.04 0.88 3.99 0.82 3.83 0.79 3.83 0.78 3.80

Molecular Magnets: Natural Orbitals versus Molecular Orbitals

Two binuclear copper diradicals, CUAQAC02 and CITLAT (Figure S1), are included in the SF-TDDFT wave function analysis benchmark study. Their experimental exchange-coupling constants (equal to ∆E within the Heisenberg-Dirac-van-Vleck model, which assumes weak interactions between the unpaired electrons) are -286 and 113 cm−1 , respectively.47, 48 This example illustrates that a compact representation of the wave function in the basis of natural orbitals provides a more facile interpretation of the underlying electronic structure than the canonical MOs. Figure 6 shows singly occupied natural and canonical MOs for the triplet states of the two binuclear copper diradicals alongside their spin-difference densities. The highest canonical MOs in CUAQAC02 and CITLAT (which have 202 and 278 electrons, respectively) appear to be delocalized. Furthermore, despite low spin-contamination of the Kohn-Sham triplet reference, α and β MOs cannot be easily matched. Thus, by considering canonical MOs only, it is difficult to ascribe overall orbital character or localization to 12


Page 12 of 25

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A

CH2

D

!"s = 1.79 (0.01) !"t = 0.99 (0.00)

!"s = 0.14 (0.01) !"t = 0.99 (0.02)

!"s = 0.19 (0.01) !"t = 0.99 (0.03)

!"s = 0.60 (0.04) !"t = 0.98 (0.1)

!"s = 0.99 (0.04) !"t = 1.00 (0.11)

!"s = 0.19 (0.00) !"t = 0.99 (0.02)

!"s = 1.85 (0.02) !"t = 1.00 (0.01)

!"s = 1.80 (0.02) !"t = 1.00 (0.02)

!"s = 1.38 (0.02) !"t = 1.00 (0.03)

!"s = 1.00 (0.09) !"t = 1.00 (0.04)

F Doublet

C

B

H

G Quartet

Doublet

! $ = 0.99 (0.46)

I

Quartet

Doublet

Quartet

Doublet

Quartet

!" = 0.99 (0.35)

! $ = 1.00 (0.49)

!" = 1.00 (0.35)

! $ = 1.00 (0.37)

! $ = 1.00 (0.40)

!" = 0.99 (0.34)

! $ = 1.00 (0.80)

!" = 1.00 (0.40)

!" = 1.00 (0.74)

! $ = 1.00 (0.36)

!" = 0.99 (0.38)

! $ = 0.99 (0.25)

!" = 1.00 (0.33)

! $ = 1.00 (0.68)

! $ = 1.00 (0.30)

!" =1.00 (0.99)

!" = 0.99 (0.99) ! $ = 0.98 (0.38) !" = 1.85 (0.02)

!" = 1.80 (0.02) !" = 0.99 (0.21)

Figure 5: EOMSF-SF-CCSD/cc-pVTZ natural orbitals of lowest singlet/doublet and triplet/quartet (low-spin) states of A-I. n ¯ = |nα + nβ |, with ∆n = |nα − nβ | provided in parentheses. n ¯ s and n ¯ t correspond to n ¯ values obtained from the occupancies of the singlet and triplet natural orbitals, respectively.

the unpaired electrons that give rise to the diradical states. In contrast, the frontier natural orbitals obtained from the triplet density matrix have expected shapes that can be described as dxy or dyz orbitals localized on the two copper centers. The spin-density differences (shown in Fig. 6) are consistent with the shapes of frontier natural orbitals.

13



CUAQAC02

CITLAT

Spin Density

Spin Density

Page 14 of 25

SOMOs

SONOs

SOMOs

SONOs

!2

!2

!2

!2

!1

!1

!1

!1

Figure 6: Spin-difference densities, unrestricted singly occupied (SO) molecular and natural frontier orbitals of high-spin triplet states of CUAQAC02 (left) and CITLAT (right) computed by B5050LYP/cc-pVTZ.

14


Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.5


Comparison between EOM-SF-CCSD, SF-TDDFT, and SF-ADC(2)s wave functions

Tables 2, 3, and 5, and 4 show DFT and ADC(2) energy gaps and state properties for the low-lying spin-states of di- and triradical benzynes, methylene, CUAQAC02, and CITLAT. nu and |nα − nβ | for triradicals are provided in the SI. Figures S4-S7 in SI show natural orbitals of the relevant spin-states at each level of theory. Diradicals with states affected by spin-contamination typically have large |nα − nβ | values. There is agreement in natural orbital character—and sometimes energetics—among DFT methods, and reasonable agreement with EOM-SF-CCSD for organic radicals. LDA and B97 fail for binuclear copper radicals with respect to energies, overestimating the degree of closed-shell character of the low-spin states relative to PBE50 and B5050LYP. The two non-hybrid functionals appear sufficient for describing the underlying wave function character of di- and triradicals with respect to the shapes of natural orbital, with the notable exception of CUAQAC02 and the B97 functional, where the dyz orbitals predicted by EOM-CCSD, PBE50, B5050LYP, and LDA are not the frontier natural orbitals (B97 frontier natural orbitals of CUAQAC02 and CITLAT shown in Figure S4 in SI). LDA and B97 appear to systematically underestimate diradical character of singlet states as reflected by lower nu,nl and n ¯ s values, relative to PBE50 and B5050LYP. SF-ADC(2) has been used to compute the energy splittings and Head-Gordon’s indices of CH2 and A-D (Tables 2 and 5). For these molecules, the singlet-triplet energy splittings (∆E) computed with SF-ADC(2) agree very favorably with the ones obtained at the EOM-SF-CCSD level of theory. The largest deviation is 0.08 eV found for CH2 . Inspecting Head-Gordon’s index nu,nl for these molecules (Table 2), we note a good overall agreement with the EOM-SF-CCSD values. The SF-ADC(2) values are consistently slightly larger by about 0.15. Further detailed investigations of the quality of the wave functions obtained with SF-ADC approaches, including higher orders of perturbation theory, will be the topic of future work. Natural orbitals of the binuclear copper complexes are of dxy and dyz type and well-localized on copper centers. The unpaired electrons of CITLAT (the copper diradical with a positive ∆E) exhibit some through-space anti-bonding interaction with p-orbitals of neighboring oxygen atoms in the bonding plane. This interaction is absent in CUAQAC02 (the copper diradical with a negative ∆E), in which singly occupied natural orbitals can be described as the dyz orbitals that give rise to δ-δ ∗ type interactions and bond orders greater than three between transition metal nuclei in more strongly interacting bimetallic complexes.59

3.6

Adiabatic Singlet-Triplet and Doublet-Quartet Gaps

All energy differences between low-lying spin states reported here are vertical energy gaps, computed at the equilibrium geometries of the ground state by optimizing the lowest energy SF-TDDFT/cc-pVTZ state. In order to assess the absolute accuracy of the methods by comparison against experimental values, one needs to consider adiabatic gaps, as was done in previous benchmark studies.28, 30, 31, 60 As a guidance for readers, here we provide a compilation of adiabatic energy gaps for the benchmark compounds considered here (with the exception of D, CUAQAC02, and CITLAT) reported elsewhere28, 30, 31, 60 and are summarized in Tables S5 and S6 in SI. Note that in the reported experimental values zero-point energies are subtracted. As one can see, EOM-SF-CCSD(dT) is within 1 kcal/mol from the experiment. EOM-SF-CCSD 15



Table 2: SF-TDDFT and ADC(2)-s energy splittings (eV) and wave function properties of lowest singlet and Ms = 0 triplet states of diradicals. Molecule

CH2

A

B

C

D

Method

∆E

PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2)

1.03 0.30 0.93 0.23 1.02 -2.38 -2.66 -2.86 -2.84 -2.46 -1.85 -2.04 -2.23 -2.22 -1.97 -0.19 -0.23 -0.52 -0.40 -0.24 0.24 0.16 0.20 0.44 0.20

nu,nl Singlet Triplet 0.30 2.00 0.08 2.00 0.33 2.00 0.04 2.00 0.33 2.01 0.09 2.00 0.04 2.00 0.01 2.00 0.01 2.00 0.32 2.17 0.19 2.00 0.12 2.00 0.01 2.00 0.03 2.00 0.41 2.17 1.47 2.00 1.29 2.00 0.30 2.00 0.63 2.00 1.58 2.18 2.02 2.02 2.01 2.02 2.00 2.00 2.00 2.00 2.24 2.24

hS 2 i Singlet Triplet 0.01 2.02 0.01 2.01 0.01 2.01 0.01 2.00 – – 0.03 2.02 0.02 2.02 0.01 2.02 0.02 2.01 – – 0.07 2.11 0.05 2.07 0.02 2.03 0.03 2.05 – – 0.03 2.03 0.02 2.02 0.01 2.02 0.02 2.02 – – 0.33 2.10 0.23 2.10 0.02 2.03 0.66 1.47 – –

delivers consistent performance and is very close to EOM-SF-CCSD(dT). This comparison justifies using EOM-SF-CCSD vertical gaps as the benchmark for SF-DFT and ADC(2).

16


Page 16 of 25

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: SF-TDDFT/cc-pVTZ energy splittings (eV) and wave function properties of lowest Ms = 1/2 doublet and quartet states of triradicals. Results for four density functionals are shown. Molecule

F

G

H

I

Method

∆E

PBE50 B5050LYP LDA B97 PBE50 B5050LYP LDA B97 PBE50 B5050LYP LDA B97 PBE50 B5050LYP LDA B97

- 2.18 -2.44 -2.80 -2.83 -2.64 -2.82 -3.16 -3.16 0.39 0.23 0.28 0.46 -0.38 -0.14 0.03 0.19

nu,nl Doublet Quartet 1.17 3.01 1.09 3.00 1.01 3.00 1.02 3.00 1.09 3.00 1.06 2.01 1.01 3.00 1.01 3.00 3.06 3.05 3.04 3.04 3.00 3.00 3.01 3.01 3.01 3.04 3.01 3.03 3.00 3.00 3.00 3.00

hS 2 i Doublet Quartet 0.87 3.96 0.84 3.88 0.78 3.80 0.80 3.81 0.79 3.80 0.79 3.79 0.76 3.75 0.78 3.77 1.42 3.81 0.99 4.10 0.79 3.80 1.12 3.64 1.79 3.21 1.77 3.16 1.25 3.29 1.74 2.91

Table 4: SF-TDDFT energy splittings (cm−1 ) and wave function properties of lowest singlet and Ms = 0 triplet states of binuclear copper diradicals. nu,nl hS 2 i Molecule Method ∆E Singlet Triplet Singlet Triplet PBE50 -123 1.97 2.00 0.01 2.01 B5050LYP -148 1.96 2.00 0.01 2.01 CUAQAC02 LDA -1420 0.47 2.00 0.01 2.01 B97 71 2.59 2.69 1.07 1.77 PBE50 100 2.00 2.00 0.01 2.02 B5050LYP 77 2.00 2.00 0.01 2.01 CITLAT LDA 411 1.96 2.00 0.01 2.01 B97 228 1.98 2.00 0.01 2.01

17



Page 18 of 25

Table 5: nu and nu,nl at the SF-TDDFT and ADC(2)-s/cc-pVTZ levels of theory. Values are provided for the lowest singlet and Ms = 0 triplet states of each diradical benchmark system. Results for four density functionals are compared. Molecule

CH2

A

B

C

D

CUAQAC02

CITLAT

Method PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 ADC(2) PBE50 B5050LYP LDA B97 PBE50 B5050LYP LDA B97

nu 2.014 2.010 2.008 2.013 2.114 2.016 2.012 2.010 2.014 3.133 2.070 2.045 2.020 2.033 3.146 2.019 2.014 2.010 2.014 3.159 2.161 2.142 2.019 2.055 3.796 2.010 2.009 2.004 2.815 2.011 2.010 2.006 2.009

Triplet nu,nl |nu − nu,nl | 2.000 0.014 2.000 0.010 2.000 0.008 2.000 0.013 2.009 0.105 2.000 0.016 2.000 0.012 2.000 0.010 2.000 0.014 2.171 0.962 2.003 0.067 2.001 0.044 2.000 0.020 2.001 0.032 2.174 0.972 2.000 0.019 2.000 0.014 2.000 0.010 2.000 0.014 2.179 0.980 2.019 0.142 2.015 0.127 2.000 0.019 2.002 0.053 2.244 1.552 2.000 0.010 2.000 0.009 2.000 0.004 2.685 0.130 2.000 0.011 2.000 0.010 2.000 0.006 2.000 0.009

18


nu 0.447 0.217 0.469 0.166 0.578 0.246 0.171 0.070 0.095 1.458 0.399 0.302 0.092 0.154 1.546 1.272 1.129 0.450 0.693 2.377 2.164 2.118 2.018 2.047 3.785 1.845 1.814 0.567 2.754 1.964 1.951 1.804 1.861

Singlet nu,nl |nu − nu,nl | 0.303 0.144 0.079 0.138 0.332 0.137 0.045 0.121 0.326 0.252 0.091 0.155 0.044 0.127 0.007 0.063 0.013 0.082 0.323 1.135 0.192 0.207 0.116 0.186 0.011 0.081 0.030 0.124 0.406 1.140 1.470 0.198 1.286 0.157 0.304 0.146 0.627 0.066 1.577 0.800 2.019 0.145 2.010 0.108 2.000 0.018 2.002 0.045 2.237 1.548 1.973 0.128 1.961 0.147 0.466 0.101 2.589 0.165 1.998 0.034 1.996 0.045 1.959 0.155 1.978 0.117

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4


Conclusion

Analysis of the natural orbitals derived from the one-particle density matrices provides insight into the extent and type of radical character for a variety of di- and triradical species. This analysis affords quantitative comparisons of electronic properties beyond energy differences computed by different methods. In particular, using these tools one can compare the performance of wave function and DFT methods. In agreement with earlier benchmark studies, we observe good agreement between EOM-SF-CCSD and SF-DFT (when using recommended functionals). The comparison of natural orbitals and their occupations computed for the lowest singlet and triplet and doublet and quartet states indicates good agreement between EOM-SFCCSD and TDDFT in most cases. We observe very good agreement among DFT functionals with regard to the character of frontier natural orbitals, however, the respective occupations (and, consequently, the effective number of unpaired electrons) vary. SF-ADC(2) results agree favorably with SF-EOM-CCSD both with respect to the relative energies of singlet and triplet states as well as with respect to the corresponding wave function character of open-shell singlets, indicating its potential as benchmark method for larger molecular systems for which SF-EOM-CCSD is no longer feasible. Our study represents the first investigation of di- and triradicals focusing on understanding state characters, rather than energies alone. For systems like the large copper-containing diradicals considered here, the canonical Kohn-Sham or Hartree-Fock orbitals fail to represent the correct bonding pattern even for high-spin states. In contrast, natural frontier orbitals, their occupations, and Head-Gordon’s index allow one to obtain a clear picture of the underlying electronic structure.

Acknowledgment This work is supported by the Department of Energy through the DE-FG02-05ER15685 grant (A.I.K.).

Associated Content Supporting Information is available free of charge on the ACS Publications website at DOI: xxx. Experimental structures of CUAQAC02 and CITLAT; additional structural data for D and H; EOM-SF-CCSD energies, properties, and natural orbitals at different geometries of B; tabulated EOM-SF-CCSD and B5050LYP energies and state properties of H2 for the triplet (high-spin and low-spin) and singlet states; expanded comparison of two Head-Gordon indices for di- and triradicals; additional NO figures for various levels of theory; compilation of results for adiabatic energy gaps; B97 NOs for CUAQAC02 and CITLAT; and all calculated and/or experimental structures used in this work.

Conflicts of Interest The authors declare the following competing financial interest(s): A.I.K. is a member of the Board of Directors and a part-owner of Q-Chem, Inc. 19



†

Corresponding author: Anna Krylov, [email protected]

20


Page 20 of 25

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References [1] Salem, L.; Rowland, C. The electronic properties of diradicals Angew. Chem. Int. Ed. 1972, 11, 92–111. [2] Bonaˇcić-Kouteck´ y, V.; Kouteck´ y, J.; Michl, J. Neutral and charged biradicals, zwitterions, funnels in S1 , and proton translocation: Their role in photochemistry, photophysics, and vision Angew. Chem. Int. Ed. 1987, 26, 170–189. [3] Krylov, A.I. Triradicals J. Phys. Chem. A 2005, 109, 10638–10645. [4] Krylov, A. I. In Reviews in Comp. Chem.; Parrill, A.L., Lipkowitz, K.B., Eds., Vol. 30; J. Wiley & Sons, 2017. [5] Nanda, K. D.; Krylov, A. I. Effect of the diradical character on static polarizabilities and two-photon absorption cross-sections: A closer look with spin-flip equation-of-motion coupled-cluster singles and doubles method J. Chem. Phys. 2017, 146, 224103. [6] Weinhold, F.; Landis, C. R. Natural bond orbitals and extensions of localized bonding concepts Chem. Ed.: Res. & Pract. Eur. 2001, 2, 91–104. [7] Weinhold, F.; Landis, C.R. Discovering Chemistry With Natural Bond Orbitals; New Jersey: John Wiley & Sons, 2012. [8] Takatsuka, K.; Fueno, T.; Yamaguchi, K. Distribution of odd electrons in ground-state molecules Theor. Chim. Acta 1978, 48, 175–183. [9] Head-Gordon, M. Characterizing unpaired electrons from the one-particle density matrix Chem. Phys. Lett. 2003, 372, 508–511. [10] Luzanov, A.V. In Practical Aspects of Computational Chemistry IV; Leszczynski, J., Shukla, M., Eds.; Springer-Verlag, New York, 2016. [11] Luzanov, A.V.; Zhikol, O.A. In Practical aspects of computational chemistry I: An overview of the last two decades and current trends; Leszczynski, J., Shukla, M.K., Eds.; Springer, 2012; p. 415. [12] Plasser, F.; Wormit, M.; Dreuw, A. New tools for the systematic analysis and visualization of electronic excitations. I. Formalism J. Chem. Phys. 2014, 141, 024106–13. [13] Plasser, F.; Bäppler, S.A.; Wormit, M.; Dreuw, A. New tools for the systematic analysis and visualization of electronic excitations. II. Applications J. Chem. Phys. 2014, 141, 024107–12. [14] Luzanov, A. V.; Casanova, D.; Feng, X.; Krylov, A. I. Quantifying charge resonance and multiexciton character in coupled chromophores by charge and spin cumulant analysis J. Chem. Phys. 2015, 142, 224104. [15] Löwdin, P.-O. Quantum theory of many-particle systems. I. Physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction Phys. Rev. 1955, 97, 1474. 21



[16] Head-Gordon, M.; Grana, A. M.; Maurice, D.; White, C. A. Analysis of electronic transitions as the difference of electron attachment and detachment densities J. Phys. Chem. 1995, 99, 14261 – 14270. [17] Martin, R.L. Natural transition orbitals J. Phys. Chem. A 2003, 118, 4775–4777. [18] Luzanov, A.V.; Sukhorukov, A.A.; Umanskii, V.E. Application of transition density matrix for analysis of excited states Theor. Exp. Chem. 1976, 10, 354–361. [19] Luzanov, A.V.; Prezhdo, O.V. Analysis of multiconfigurational wave functions in terms of hole-particle distributions J. Chem. Phys. 2006, 124, 224109. [20] Surján, P. R. Natural orbitals in CIS and singular-value decomposition Chem. Phys. Lett. 2007, 439, 393–394. [21] Mayer, I. Using singular value decomposition for a compact presentation and improved interpretation of the CIS wave functions Chem. Phys. Lett. 2007, 437, 284–286. [22] Wagner, P. J.; Hammond, G. S. Mechanisms of photochemical reactions in solution. xxxviii.’ quenching of the type i1 photoelimination reaction J. Am. Chem. Soc. 1966, 88, 1245–1251. [23] Seinfeld, J. H.; Pandis, S. N. Atmospheric chemistry and physics: from air pollution to climate change; John Wiley & Sons, 2016. [24] Kahn, O. Magnetism: a supramolecular function, Vol. 484; Springer Science & Business Media, 2013. [25] Krylov, A. I. Size-consistent wave functions for bond-breaking: The equation-of-motion spin-flip model Chem. Phys. Lett. 2001, 338, 375–384. [26] Levchenko, S. V.; Krylov, A. I. Equation-of-motion spin-flip coupled-cluster model with single and double substitutions: Theory and application to cyclobutadiene J. Chem. Phys. 2004, 120, 175–185. [27] Shao, Y.; Head-Gordon, M.; Krylov, A. I. The spin-flip approach within time-dependent density functional theory: Theory and applications to diradicals J. Chem. Phys. 2003, 118, 4807–4818. [28] Bernard, Y. A.; Shao, Y.; Krylov, A. I. General formulation of spin-flip time-dependent density functional theory using non-collinear kernels: Theory, implementation, and benchmarks J. Chem. Phys. 2012, 136, 204103. [29] Lefrancois, D.; Wormit, M.; Dreuw, A. Adapting algebraic diagrammatic construction schemes for the polarization propagator to problems with multi-reference electronic ground states exploiting the spin-flip ansatz J. Chem. Phys. 2015, 143, 124107. [30] Slipchenko, L. V.; Krylov, A. I. Singlet-triplet gaps in diradicals by the spin-flip approach: A benchmark study J. Chem. Phys. 2002, 117, 4694–4708.

22


Page 22 of 25

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[31] Slipchenko, L. V.; Krylov, A. I. Electronic structure of the trimethylenemethane diradical in its ground and electronically excited states: Bonding, equilibrium structures and vibrational frequencies J. Chem. Phys. 2003, 118, 6874–6883. [32] Slipchenko, L. V.; Krylov, A. I. Electronic structure of the 1,3,5,-tridehydrobenzene triradical in its ground and excited states J. Chem. Phys. 2003, 118, 9614–9622. [33] Slipchenko, L. V.; Munsch, T. E.; Wenthold, P. G.; Krylov, A. I. 5-dehydro-1,3quinodimethane: A hydrocarbon with an open-shell doublet ground state Angew. Chem. Int. Ed. 2004, 43, 742–745. [34] Vanovschi, V.; Krylov, A. I.; Wenthold, P. G. Structure, vibrational frequencies, ionization energies, and photoelectron spectrum of the para-benzyne radical anion Theor. Chim. Acta 2008, 120, 45–58. [35] Hossain, E.; Deng, S. M.; Gozem, S.; Krylov, A. I.; Wang, X.-B.; Wenthold, P. G. Photoelectron spectroscopy study of quinonimides J. Am. Chem. Soc. 2017, 139, 11138–11148. [36] Lefrancois, D.; Rehn, D.; Dreuw, A. Accurate adiabatic singlet-triplet gaps in atoms and molecules employing the third- order spin-flip algebraic diagrammatic construction scheme for the polarization propagator J. Chem. Phys. 2016, 145, 084102. [37] Lefrancois, D.; Tuna, D.; Martinez, T.J.; Dreuw, A. The spin-flip variant of the algebraicdiagrammatic construction yields the correct topology of S1 /S0 conical intersections J. Chem. Theory Comput. 2017, 13, 4436–4441. [38] Manohar, P. U.; Krylov, A. I. A non-iterative perturbative triples correction for the spinflipping and spin-conserving equation-of-motion coupled-cluster methods with single and double substitutions J. Chem. Phys. 2008, 129, 194105. [39] Luxon, A.; Orms, N.; Kanters, R.; Krylov, A. I.; Parish, C. An ab initio exploration of the bergman cyclization J. Phys. Chem. A 2017; inpress. [40] Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular electronic structure theory; Wiley & Sons, 2000. [41] Krylov, A. I. The spin-flip equation-of-motion coupled-cluster electronic structure method for a description of excited states, bond-breaking, diradicals, and triradicals Acc. Chem. Res. 2006, 39, 83–91. [42] Sears, J. S.; Sherrill, C. D.; Krylov, A. I. A spin-complete version of the spin-flip approach to bond breaking: What is the impact of obtaining spin eigenfunctions? J. Chem. Phys. 2003, 118, 9084–9094. [43] Dreuw, A.; Wormit, M. The algebraic diagrammatic construction scheme for the polarization propagator for the calculation of excited states WIREs Comput. Mol. Sci. 2015, 5, 82–95. [44] Krylov, A. I.; Sherrill, C. D. Perturbative corrections to the equation-of-motion spinflip SCF model: Application to bond-breaking and equilibrium properties of diradicals J. Chem. Phys. 2002, 116, 3194–3203. 23



[45] Christiansen, O.; Koch, H.; Jørgensen, P. The 2nd-order approximate coupled-cluster singles and doubles model CC2 Chem. Phys. Lett. 1995, 243, 409–418. [46] Sherrill, C.D.; Leininger, M.L.; Van Huis, T.J.; Schaefer III, H.F. Structures and vibrational frequencies in the full configuration interaction limit: Predictions for four electronic states of methylene using triple-zeta plus double polarization (TZ2P) basis J. Chem. Phys. 1998, 108, 1040–1049. [47] deMeester, P.; Fletcher, S. R.; Skapski, A. C. Refined crystal structure of tetra-p-acetatobisaquodicopper(l1) J. C. S. Dalton 1973, 23, 2575–2578. [48] Youngme, S.; Phatchimkun, J.; Wannarit, N.; Chaichit, N.; Meejoo, S.; vanAlbada, G.A.; Reedijk, J. Refined crystal structure of tetra-p-acetato-bisaquodicopper(l1) Polyhedron 2008, 27, 304–318. [49] Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple Phys. Rev. Lett. 1996, 77, 3865–3868. [50] Becke, A. D. A multicenter numerical integration scheme for polyatomic molecules J. Chem. Phys. 1988, 88, 2547. [51] Lee, C.; Yang, W.; Parr, R.G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron-density Phys. Rev. B 1988, 37, 785–789. [52] Becke, A.D. Density-functional thermochemistry. V. Systematic optimization of exchangecorrelation functionals J. Chem. Phys. 1997, 107, 8554–8560. [53] Wang, F.; Ziegler, T. Time-dependent density functional theory based on a noncollinear formulation of the exchange-correlation poential J. Chem. Phys. 2004, 121, 12191. [54] Vahtras, O.; Rinkevicius, Z. General excitations in time-dependent density functional theory J. Chem. Phys. 2007, 126, 114101. [55] Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A.T.B.; Wormit, M.; Kussmann, J.; Lange, A.W.; Behn, A.; Deng, J.; Feng, X., et al. Advances in molecular quantum chemistry contained in the Q-Chem 4 program package Mol. Phys. 2015, 113, 184–215. [56] Krylov, A. I.; Gill, P. M. W. Q-Chem: An engine for innovation WIREs Comput. Mol. Sci. 2013, 3, 317–326. [57] Cristian, A. M. C.; Shao, Y.; Krylov, A. I. Bonding patterns in benzene triradicals from structural, spectroscopic, and thermochemical perspectives J. Phys. Chem. A 2004, 108, 6581–6588. [58] Smith, C.E.; Crawford, T.D.; Cremer, D. The structures of m-benzyne and tetrafluorom-benzyne J. Chem. Phys. 2005, 122, 174309. [59] Wagner, F. R.; Noor, A.; Kempe, R. Ultrashort metalmetal distances and extreme bond orders Nat. Chem. 2009, 1, 529–536. [60] Wang, T.; Krylov, A. I. Electronic structure of the two dehydro-meta-xylylene triradicals and their derivatives Chem. Phys. Lett. 2006, 426, 196–200. 24


Page 24 of 25

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC graphics

𝑛"occ = 0.14

𝑛"occ = 0.19

𝑛"occ = 0.60

𝑛"occ = 1.85

𝑛"occ = 1.80

𝑛"occ = 1.38

25


Characterizing Bonding Patterns in Diradicals and ... - ACS Publications

Recommend Documents