Assignment of the 6,6â²-Dioxo-3,3â²-biverdazyl Ground State by

2010

J. Phys. Chem. A 2010, 114, 2010–2021

Assignment of the 6,6′-Dioxo-3,3′-biverdazyl Ground State by Using Broken Symmetry and Spectroscopy Oriented Configuration Interaction Techniques Saba M. Mattar* and Hisham M. Dokainish UniVersity of New Brunswick, Department of Chemistry and Centre for Laser, Atomic, and Molecular Sciences, Fredericton, New Brunswick, Canada E3B 6E2 ReceiVed: NoVember 8, 2009

The singlet-triplet energy differences (∆EST) and the Heisenberg-Dirac-van-Vleck exchange parameter (J) of 6,6′-dioxo-3,3′-biverdazyl (BVD) have been studied by hybrid density functional (HDF), broken symmetry (BS), and spectroscopy oriented configuration interaction (SORCI) methods. Energy surface scans as a function of the dihedral angle between the two verdazyl rings (φN2C3C3′N2′) have been performed. The BS computations predict an antiferromagnetic ground state. However, the diradical index (RBS) ranges from 97.5 to 99.9%, implying that the interactions between the two unpaired electrons are very weak. To calculate J and ∆EST, the multireference character introduced by these weak spin-spin interactions must be taken into account. Consequently, multireference difference dedicated configuration interaction (MRDDCI) methods, as implemented in the SORCI procedure, are used. The in-plane π (IPπ), out-of-plane π (OPπ), and σ configurations are included in the CI expansions in a balanced fashion. The OPπ-OPπ and OPπ-IPπ overlaps are the predominant factors that influence the J and ∆EST as a function of φN2C3C3′N2′ and cause them to peak around 40 and 140°. In these regions, the antiferromagnetic interactions are minimal, and the MRDDCI methods predict a triplet ground state. At φN2C3C3′N2′ ) 0, ∆EST[MRDDCI3(14,12)] is in excellent agreement with that of 1,1′,5,5′-tetramethyl-6,6′-dioxo-3,3′-biverdazyl determined experimentally from electron paramagnetic resonance (EPR) spectroscopy and differs only by 2.3%. Furthermore, ∆EST[MRDDCI3(14,12)] is consistently smaller than JY as the verdazyl rings rotate with respect to each other. This corroborates the theory that the HDF-BS technique increases the singlet-triplet energy gap and favors the singlet state. Because the SORCI method is specifically designed for large molecules, the present very good results open the way for the computation of the magnetic properties of larger molecules by the SORCI method. To the best of our knowledge, this is the first time that ∆EST has been computed by the MRDDCI3 method by utilizing such a large CI reference space for a molecule of this size. 1. Introduction The coupling of unpaired valence electrons in organic solids plays a crucial role in determining their magnetic, semiconducting, conducting, and superconducting properties.1-4 The simplest organic solids that undergo such intramolecular electronic interactions are diradicals that possess two unpaired electrons. The two electrons may be localized in specific regions of the molecule or delocalized over large portions.5 The study of these diradicals may help understand the theoretical principles that govern the spin multiplicity (singlet or triplet) of their ground and low-lying excited states. These, in turn, would facilitate the correlation between their structures and their bonding and molecular properties.6 The most intuitive way to form a diradical is to bond two radicals to one another such that their unpaired electrons remain separated and spatially confined to their original moieties. This may be accomplished by connecting the two radicals via two atoms that have a node in their singly occupied molecular orbitals (SOMOs). The resulting two molecular orbitals, formed from the two radical SOMOs, should be nonbonding with respect to the linking atoms and sever the link between the two unpaired electrons. However, if additional interactions between the two unpaired electrons should occur, one should not assume * Corresponding author. Tel.: +1(506) 447 3091; Fax: +1(506) 453 4981; E-mail: [email protected].

that these interactions are necessarily antiferromagnetic and result in a singlet ground state. For example, tetramethylene ethane7 with D2d symmetry, the slightly twisted near planar C2 diradical 2,3-bis(methylene)-1,3-cyclohexadiene,8 and 9,9′bianthryl9 3 all have a triplet ground states. To further complicate matters, bis(nitronyl nitroxide), which is in some respects an analogue of 9,9′-bianthryl, is found to have a singlet ground state.10,11 Thus, the factors that determine whether a molecule containing two unpaired electrons is a diradical, singlet, or triplet state are still unclear and need further investigation. Bis-verdazyls are a class of molecules that fall in the above category. In general, verdazyls are monoradicals or can be linked together to form di- and triradicals and so forth. The monoverdazyl base units are π-radicals, as shown in Scheme 1. Therefore, poly verdazyls are potential building blocks for conducting devices or molecular magnets.12 A large number of verdazyl radicals have been synthesized and studied by electron paramagnetic resonance (EPR) spectroscopy. For example, among the most recent ones is the tetrathiafulvalene (TTF) radical.13 To synthesize verdazyl radicals that are efficient molecular magnets or spintronic devices, it is imperative to comprehend their self-assembly and aggregation mechanisms in the solid state.12 It is also just as important to first understand how their electronic structure and bonding affect their individual magnetic

10.1021/jp910643e  2010 American Chemical Society Published on Web 01/05/2010

6,6′-Dioxo-3,3′-biverdazyl Ground State with BS and SORCI Techniques

J. Phys. Chem. A, Vol. 114, No. 4, 2010 2011

SCHEME 1

and electronic properties. This is accomplished by relating their experimental spin Hamiltonian parameters, such as g,14-16 fine tensor (D), and nuclear hyperfine (A) tensor components17 and the Heisenberg-Dirac-van-Vleck (HDvV) coupling constant (J) with the corresponding computed values. A good correlation between theory and experiment gives a clear picture of the radical’s structure and bonding together with its total and net spin density distributions. Verdazyl monoradicals are characterized by ground-state wave functions that are predominantly single Slater Kohn-Sham determinants. As a result, hybrid density functional (HDF) methods, combined with moderate basis sets, should be very successful in calculating their electronic structure, optimal geometries, and A and g tensors.18 The next step is to try and understand the structure and bonding properties of verdazyl diradicals. In particular, one needs to determine the strength and nature of the interaction between the two unpaired electrons. A ferromagnetic interaction leads to a triplet state, whereas an antiferromagnetic one results in either a closed-shell singlet or an open-shell diradical. The 1,1′,5,5′-tetramethyl-6,6′-dioxo-3,3′-biverdazyl diradical (TMBVD) has an unpaired electron on each ring. Electron paramagnetic resonance studies have shown that these two electrons couple antiferromagentically, leading to a singlet ground state. In solution, its singlet ground state lies 758 cm-1 below the triplet state, whereas in the solid state, it is 887 cm-1.6 The TMBVD crystal structure indicates that in the solid state, it has local D2h symmetry in the vicinity of its C3sC3′ linking bond. However, in solution, it is twisted with local D2 symmetry.6 In general, the theoretical prediction whether a molecule has a singlet or triplet ground state is not an easy task, particularly for large molecules. Accurate calculation of the singlet-triplet energy difference, ∆EST, is even more difficult and optimally requires extensive multireference-configuration interaction (MRCI) calculations.19 The 6,6′-dioxo-3,3′-biverdazyl (BVD), shown in Figure 1, is an analogue of TMBVD where the four methyl groups are replaced by hydrogens. It is a benchmark molecule, and its ∆EST has been previously calculated by using various computational methods. Chung and Lee20 used a multiconfiguration selfconsistent field (MCSCF) method as well as the B3LYP variant of the HDF technique. Barone et al.18 estimated ∆EST and J of BVD by using HDF theory combined with broken-symmetry (BS) solutions of its singlet state. To take into consideration solvent effects, the polarizable continuum model (PCM) was also used.18 Recently, the complete-active-space second-order perturbation (CASPT2) method has also been used in an attempt to get a more accurate estimation of ∆EST.21,22 In this article, the HDF-BS, complete-active space selfconsistent field (CASSCF), and multireference-difference-

Figure 1. Structure and atomic numbering of the 6,6′-dioxo-3,3′biverdazyl (BVD) in its planar D2h and orthogonal D2d conformations.

dedicated-configuration-interaction (MRDDCI) techniques,23 in the form of Neese’s spectroscopy-oriented-configuration-interaction (SORCI) procedure,24,25 are used to calculate ∆EST and J as a function of the N2C3C3′N2′ dihedral angle (φN2C3C3′N2′) from 0 to 180°. The paper is partitioned as follows. In Section 2, a brief theoretical background of the required techniques is presented. This is followed in Section 3 by the computational details. The results and discussion, in Section 4, are broken down into different subsections. The first analyzes the resulting HDF-BS wave functions and their diradical indices. The next subsection gives a concise explanation why the introduction of static correlation in the CASSCF procedure reduces the magnitudes of ∆EST and J compared to those obtained by the HDF-BS method. The next subsection describes the introduction of dynamic correlation, in a systematic fashion. This includes the MRDDCI methods23 as incorporated in the SORCI procedure.24 The behavior of ∆EST as a function of φN2C3C3′N2′ is then interpreted in terms of ferro- and antiferromagnetic spin interaction pathways in conjunction with spin polarization. Finally, a summary of the results and conclusions are given. 2. Theoretical Background The magnetic exchange interaction between two spin systems is given by the HDvV Hamiltonian

HHDvV ) -2JSASB

(1)

where SA and SB are the spin angular momentum operators and J is the exchange coupling constant. In the BVD case, SA ) 1/2 and SB ) 1/2, which leads to a triply degenerate high-spin (HS) state with SMax ) SA + SB ) 1 and a low-spin (LS) state where SMin ) SA - SB ) 0. A positive sign of J indicates a ferromagnetic ground state, whereas a negative sign indicates an antiferromagnetic one.26,27 The single determinant BS method is one of the most popular techniques used to estimate J for large molecules because of its speed and simplicity.28 It depends on evaluating J by determining the triplet HS and the singlet BS total energies. The HS state energy can be determined by an unrestricted SCF procedure. On the other hand, the BS wave functions are single determinants, similar to |vV| or |Vv|, and are not eigenfunctions of the total spin operator. Consequently, the BS solution cannot describe the true ground singlet state which requires a linear combination of at least two determinants, |vV| and |Vv|. Nevertheless, mathematical relations exist between the HS and BS

2012

J. Phys. Chem. A, Vol. 114, No. 4, 2010

Mattar and Dokainish

energies and wave functions that enable one to estimate J. Two relations are normally used. The first is due to Noodleman28,29 where J, for the weak coupling case, is

JN ) -

EHS - EBS 2 SMax

(2)

Here, EHS and EBS represent the energy eigenvalues of the HS and BS states, respectively. A second more general equation was later suggested by Yamaguchi who derived the expression for J in terms of the spin expectation values30

JY ) -

EHS - EBS 2 2 〈SˆMax 〉 - 〈SˆBS 〉

(3)

Although the HDF-BS method is simple and gives adequate J values in some diradical systems, it usually overestimates the stability of the singlet state, leading to negative J values that are too large. Modifications have been proposed to improve the original HDF-BS formula, such as the constrained HDF-BS approach proposed by the van Voorhis group.31 It gives better results compared to the conventional HDF-BS methods for organometallic and inorganic metal clusters. However, it is used only if the spins are well localized and weakly coupled.31 To obtain accurate ∆EST and J values by MRCI methods, good-quality initial wave functions are needed. These wave functions should preferably include static correlation. They can be obtained from CASSCF calculations.32-34 The CASSCF wave function, |ΨIs〉, for the Ith state of spin S consists of a linear combination of configuration state functions (CSFs)

|Ψls〉 )

∑ Ckl|Φsk〉

(4)

k

The initial reference space in the BVD case is constructed from the CASSCF wave functions, |ΨIs〉. The classification of orbitals as frozen, internal, active, or external is crucial. Frozen orbitals are totally excluded from the CI calculations. Internal orbitals are always left doubly occupied when constructing the reference space. The active orbitals accommodate the active electrons and are partially occupied in at least one reference. The rest of the orbitals are considered as external and are empty in all references. If n is the number of active electrons and m is the number of active orbitals, the resulting complete reference space is abbreviated as CAS(n,m). Throughout this article, internal orbitals are labeled i, j, and so forth, active orbitals are, p, q, r, s, and so forth, and external orbitals are a, b, and so forth.24 The active reference space is built up from the orbitals around the HOMO-LUMO gap. Therefore, it is important that these orbitals possess the right type of bonding characteristics. This is accomplished by examining these orbitals and properly ordering them before starting the SORCI calculations. The computational times and sizes of the calculations are very sensitive to the number of configurations in the reference space. To make these calculations practical, the reference space should be as small as possible. Initially, the Hamiltonian is diagonalized in a new basis of CSFs formed by spin and space adapting the original orbitals. If the weights of the resulting roots are less than a specific threshold, Tpre, they are discarded.24 The new selected and reduced reference space, St{Φ}, is then used to construct the multiconfigurational zero-order wave function,

∑ cνiφi(r)

(5)

(6)

µεS

for all the required I states. The diagonalization of the Hamiltonian with respect to this space yields the zeroth-order energies EI(0). A Mo¨ller-Plesset multireference perturbation (MR-MP2) computation, for each Ith state, is carried, out and the secondorder energy is obtained from

Each |Φks〉 is formed from a set of orthonormal molecular orbitals, ψi(r), that are spin and space symmetry adapted to the desired spin state, S. In turn, every individual ψi(r) is expanded in terms of a set of basis functions, φν(r) as

ψi(r) )

∑ cµI(0)Φµ

Ψ(0) I )

E(2) I

)

2 ˆ |〈Ψ(0) I |H |Φµ〉|

∑ 〈Ψ(0)|Hˆ |Ψ(0)〉 - 〈Φ |Hˆ |Φ 〉

µ∉S

I

0

I

µ

0

(7)

µ

where the zeroth-order Hamiltonian, obtained by using the second quantization formalism, has the explicit form24

ν

The CASSCF energy is variationally minimized by adjusting the cνi and Ckl coefficients. At convergence, the energy gradients with respect to cνi and Ckl are all approximately zero. Although the CASSCF method is an enhancement over conventional SCF techniques, it lacks dynamic correlation and is still insufficient to accurately calculate ∆EST and J. Dynamic correlation is added via multireference methods. These kinds of calculation are very expensive and time-consuming. A practical alternative is to use truncated MRCI expansions. Here, Neese’s SORCI multistep procedure is an appropriate MRCI method designed to efficiently calculate energy differences between electronic spin states.24 The ORCA computer program, that implements the SORCI procedure, is extremely flexible and allows the user to have complete control over the calculations in a systematic way. A detailed description of the method and choices is given in a flowchart in the original article and manual.24,35 Here, only the relevant theory will be concisely described.

ˆ0 ) H

∑ {〈ψp|hˆ|ψp〉 + ∑ 〈Ψ(0)I | ∑ arσ+ asσ|Ψ(0)I 〉 × p

[∫

r,s

σ

ψp(r1) ψq(r1) ψr(r2) ψs(r2) dr1 dr2 |r1 - r2 | 1 ψp(r1) ψr(r1) ψq(r2) ψs(r2) dr1 dr2 + ap ap 2 |r1 - r2 |

∫

]

(8)

Here, hˆ is the one-electron Hamiltonian that includes all frozen and core MOs. The a+ and a are the Fermion creation and annihilation operators, respectively. The choice of the states to be included in eq 7 is important. Configurations resulting from the excitations of two electrons in the internal space to two orbitals in the external space, |ψa,b i,j 〉, are not included. Their effect is minimal because the energy differences due to these configurations approximately cancel each other.24 This is the basic premise behind the difference

6,6′-Dioxo-3,3′-biverdazyl Ground State with BS and SORCI Techniques dedicated CI technique first suggested by Malrieu, Caballol, and co-workers.23 All single excitations of the type |ψpi 〉, |ψpa〉, |ψai 〉, and |ψpq〉 are included because they are important in the accurate calculation of the spin polarization. The double excitations that leave the r,a a,b internal orbitals intact, |ψr,s p,q〉, |ψp,q〉, and |ψp,q〉, and those in which 〉, are included next. Adding the external orbitals are empty, |ψp,q i,j q,r q,a 〉 and |ψi,p 〉, that involve the the configurations of the form |ψi,p internal, active, and external spaces, results in summations that grow as the square of the molecular size.23 Thus, it is named a MRDDCI2 calculation.24 To further increase the amount of dynamic correlation and accuracy, one incorporates in eq 7 the configurations that include double excitations from the internal space to the active and p,a 〉. In addition, the double excitations from external spaces, |ψi,j a,b 〉, are also the internal and active space to the external one, |ψi,j 24 included. These additions cause the expansions to increase as the third power of the molecular size and hence are classified as MRDDCI3. The aim of the MRPT2 calculations is not only to determine the second-order corrections to the energy in eq 7 but also to classify the interactions arising from the first-order interacting space {R}, as strong {R′} or weak {R″}.24 If an individual energy correction term in eq 7 is larger than a specified threshold Tsel,

Tsel e

2 ˆ |〈Ψ(0) I |H |Φµ〉|

(9)

(0) ˆ ˆ 〈Ψ(0) I |H0 |ΨI 〉 - 〈Φµ |H0 |Φµ〉

it is considered as a strong perturber, and the corresponding Φµ is placed in the R′ space. If this condition is not fulfilled, it is a weak perturber and is included in the R″ space.24 The total Hamiltonian is then diagonalized in the combined {S + R′} CSF space to yield the energy eigenvalues EI(a) and the corresponding ΨI(a) states. This procedure automatically eliminates any intruder states.24 From eq 7, the second-order energy contributions arising from the unselected R″ space are

ˆ |Φ 〉| 2 |〈Ψ(0) |H

µ ∑ 〈Ψ(0)|Hˆ |Ψ(0)I 〉 - 〈Φ ˆ |Φ 〉 |H

E(unsel) ) I

µ∈R″

I

0

I

µ

0

(10)

µ

The above diagonalization procedure results in unlinked contributions. This is remedied via a Davidson correction first implemented by Buenker and co-workers36,37

(

Ecorr ) 1I

∑ cµI2 )( ∑ cµI′ cνI′ 〈Φµ|Hˆ|Φν〉 - EaI )

µ∈S

µ,ν∈S

(11)

Because one has to compare different points on the energy φN2C3C3′N2′ curves, an energy correction,

∆DDCI3 )

∑

(2 - δij) ×

iej,aeb

〈ψiψa |ψjψb〉(2〈ψiψa |ψjψb〉 - 〈ψiψb |ψjψa〉) εi + εj - (εa + εb)

(12)

as first suggested by Caballol and co-workers, is warranted.38 Thus, the final total energy is

J. Phys. Chem. A, Vol. 114, No. 4, 2010 2013 (corr) EI ) E(a) + E(unsel) + ∆DDCI3 I + EI I

(13)

The ORCA algorithm also has the option of generating approximate average natural orbitals (AANOs). The reduced density matrices for a number of specified states are generated from their natural orbitals, averaged, and then diagonalized. This results in the {ψAANO} set of natural orbitals.24 The AANOs may be trimmed by a threshold value Tnat and reused as starting orbitals in the first step of ORCA program. Repeating the calculations with AANOs as starting orbitals gives a better description of the molecule’s electronic structure for all specified states.24 3. Computational Details All the electronic structure calculations were done with the ORCA suite of programs.35 Geometry optimizations performed by using HDF computations are only appropriate for molecular electronic states that are expressed by a dominant single Kohn-Sham determinant. Consequently, the BVD molecule was geometry optimized in its lowest triplet state. The initial bond distances and angles were taken from the TMBVD solidstate crystal structure.6 The relaxed energy scans and partial geometry optimizations that followed were performed by fixing the φN2C3C3′N2′ torsion angle at a particular value and optimizing the remaining bond distances and angles by using the unrestricted B1LYP and BP86 HDF methods.39,40 The geometries resulting from the relaxed energy scans are also known as the flexible rotor model (FRM).18 This process was repeated at other φN2C3C3′N2′ angles in the range of 0-180° in increments of 10°. The BS solutions of the LS state were then calculated by using these optimized geometries. Their energies were then used, in conjunction the triplet-state energies, to estimate ∆EBS and the Heisenberg-Dirac-van-Vleck exchange parameter, J. Barone’s EPR-II orbital basis sets for H, C, O, and N were employed.41 The resolution of the identity approximation (RI)42 was used in all calculations. The auxiliary basis sets, needed for the RI approximation, were of the split valence type developed by Ahlrichs.43 In addition to BS calculations, restricted configuration interaction (CI) methods were also used to calculate ∆EST. The first step starts with a set of natural orbitals obtained from a coupled electron-pair approximation (CEPA/3) because they are better than those of HDF or MP2.44-46 The CEPA reference configurations were generated from all possible excitations of two electrons in two orbitals. The resulting natural orbitals were further used as input for CASSCF calculations.34,47,48 Two sets of calculations, CASSCF(2,2) and CASSCF(6,4), were performed. The SORCI procedure,24 with the threshold values Tpre ) -4 10 , Tsel ) 10-6, and Tnat ) 10-4, was employed to account for the dynamic correlation effects. All excitations of two electrons in two orbitals, CAS(2,2), were used as references. The resulting AANOs were subsequently used in an MRDDCI3 treatment with a larger reference space, CAS(14,12), to give the final energy results. The values of Tpre, Tsel, and Tnat were kept the same in the MRDDCI3 calculations. In the strictest sense, the classification of the molecular orbitals as σ and π only applies to linear molecules of D∞h or C∞V symmetry. However, the BVD molecule may have D2h, D2, and D2d symmetry depending on the value of φN2C3C3′N2′. Table 1 lists the irreducible representations of the BVD orbitals and the bonding classification for the three symmetries. The bonding character can be qualitatively categorized as σ, ring IPπ, and

2014



TABLE 1: BVD Orbital Symmetries, Irreducible Representations (Irepp.),a and Bonding Character

TABLE 2: Calculated BVD Energy Differences and J Valuesa

D2hb orbital symmetry Irrep. ag b1g b2g b3g au b1u b2u b3u

x f (-x) sym sym antisym antisym antisym antisym sym sym

y f (-y) sym antisym sym antisym antisym sym antisym sym

z f (-z) sym antisym antisym sym antisym sym sym antisym

bonding e

σ OPπ OPπ* IPπ* OPπ* σ*e IPπ OPπ

D 2c

D2dd

Irrep.

Irrep.

a b1 b2 b3 a b1 b2 b3

a1 a2 e e b1 b2 e e

a

Sym and antisym stand for symmetric and antisymmetric with respect to a specific coordinate inversion, respectively. b Planar BVD, φN2C3C3′N2′ ) 0°. c Twisted BVD rings, 0 < φN2C3C3′N2′ < 90° and 90 < φN2C3C3′N2′ < 180°. d Orthogonal BVD rings, φN2C3C3′N2′ ) 90°. e for atoms lying along the majorC2 x-axis.

φN2C3C3′N2′

∆EBS

JN

JY

∆EST exp.

0 10 20 30 40 50 60 70 80 90

-750.09 -686.61 -551.64 -417.54 -307.56 -256.59 -253.15 -268.54 -290.48 -292.27

-750.09 -686.61 -551.64 -417.54 -307.56 -256.59 -253.15 -268.54 -290.48 -292.27

-762.32 -699.84 -565.47 -429.43 -315.60 -261.25 -254.59 -266.68 -285.75 -286.55

-887.6b

Units of JY, JN, ∆EBS, and ∆EST are in cm-1. value for the corresponding TMBVD.6 a

b

Experimental

ring OPπ. In addition, the orbitals and electrons may be bonding or antibonding with respect to the C3sC3′ inter-ring bond. The antibonding orbitals are labeled with an extra /. Classifying the bonding approximately as σ, IPπ, and OPπ simplifies the forthcoming discussion and interpretation of the results. Visualization of the one-electron orbitals and natural orbitals is essential to judiciously include the appropriate orbitals in the CI reference space. This was accomplished with the gOpen Mol program.49,50 4. Results and Discussion Irrespective of φN2C3C3′N2′, the BVD atoms can be grouped in pairs that are related to one another spatially. From Figure 1, the atom pairs in parentheses, (N1,N5′), (N2,N4′), (C3,C3′), ..., (O7,O7′), are all spatially equivalent. Therefore, the computed JY, JN, and ∆E values are symmetrical about φN2C3C3′N2′ ) 90°. For example, the exchange coupling constants at φN2C3C3′N2′ ) 80 and 100° are identical and so forth. For the sake of conciseness, only values in the range of φN2C3C3′N2′ ) 0-90° are listed in the tables. The most stable BVD configuration is obtained when φN2C3C3′N2′ ≈ 40°. This is a characteristic behavior of molecules made up of two planar rings that are rotated with respect to one another along their inter-ring connecting bond. A typical example is that of diphenyl.46 A simple explanation of the shapes of these curves comes from the fact that delocalization of the π-type electrons over the two rings is favored when φN2C3C3′N2′ ) 0°. On the other hand, the steric repulsion between the two verdazyl rings is minimized when φN2C3C3′N2′ ) 90°. Therefore, the most stable configuration, as previously noted by Barone and co-workers,18 is expected to occur between 0-90°. The delocalization of electrons over the two rings could arise from OPπ bonding such as the b3u one-electron molecular orbitals (MOs) or IPπ bonding due to b2u MOs. In the case of biphenyl, the steric repulsion between the two rings is postulated to arise from the four ortho CH moieties. In the BVD case, the N2, N2′, N5, and N5′ lone pairs are the main cause of the steric repulsion. 4.1. BS Calculations. The computed JN and JY values as a function of φN2C3C3′N2′ along the inter-ring C3sC3′ bond are listed in Table 2. The table also contains the corresponding values of ∆EBS and the experimentally determined ∆EST value for TMBVD. All computations, obtained from eqs 2 and 3, give negative JN and JY values, indicating that the triplet state is higher than the BS state. Therefore, the current HDF-BS calculations predict

Figure 2. BVD magnetic exchange parameters, JY, JN, and ∆EBS, as a function of the inter-ring dihedral angle, φN2C3C3′N2′. ∆EBS and JY are indistinguishable.

an antiferromagnetic interaction between the two verdazyl units, irrespective of the value of torsion angle. This is in accordance with all previous calculations performed for this molecule. Equation 2, used to determine JN, only applies in the weakcoupling case between the two spins, SA and SB. On the other hand, eq 3 gives a value for JY that is general and applies in cases where the spin-spin coupling may be either weak or strong. Table 2 shows that JY and JN are very close. Thus, the system is a diradical with relatively weak interactions between the two unpaired electrons. Figure 2 gives the computed values of JY, JN, and ∆EBS for φN2C3C3′N2′ ranging from 0-180°. The JY, JN, and ∆EBS appear to be superimposed on each another. The magnitudes of JY, JN, and ∆EBS decrease as φN2C3C3′N2′ is increased from 0 to 50°. However, they are still negative, indicating a singlet ground state. In the range of φN2C3C3′N2′ ) 50-90°, the values of these three parameters drop slightly. The shape of the curve and the magnitudes of JY and JN are also in good agreement with the results obtained by Barone et al. in their BVD BS study.18 In previous calculations, the BVD molecule where the four methyl groups were replaced by H atoms was used instead of TMBVD.18,20 The comparison of the BVD results with the experimental TMBVD results was justified because the substitution did not alter the optimized geometries significantly in vacuo.18 Here, the spin densities on the C3, C3′, N2, N2′, N4, and N4′ atoms of BVD and TMBVD were also compared. The maximum difference was only 0.001 e/a03, indicating that it is indeed justifiable to compare JY, JN, and ∆EBS of BVD with those of TMBVD.6 Table 2 indicates that the calculated BVD JY and JN values are smaller than those of TMBVD by approximately 18%.

6,6′-Dioxo-3,3′-biverdazyl Ground State with BS and SORCI Techniques 4.2. Diradical Indices and Nature of the BS Wave Functions. The BS wave function, |φBS〉, is obtained from an unrestricted SCF procedure that uses a single configuration. However, the BVD singlet state can only be properly represented by a multiconfiguration wave function. This renders the interpretation and visualization of |φBS〉 in terms of conventional Kohn-Sham orbitals difficult. Neese has devised a procedure for analyzing the BS solutions in terms of valence bond-type orbitals.51 Thus, the percent ionic, neutral, triplet, and diradical character of |φBS〉 may be obtained.51,52 In the simple unrestricted Hartree-Fock case of two electrons in two orbitals, the magnetic orbitals are defined in terms of R,β51 the R and β spatial molecular orbitals, ψR,β + and ψ-

ηAR,β ) ηBR,β )

1 R,β R,β ) (ψ+ + ψ√2

(14)

1 R,β R,β ) (ψ+ - ψ√2

(15)


Rβ j βi 〉.51 If ψRi ) φAR maximum similarity and overlap, SVar ) 〈ψRi |ψ β β j and ψi ) φ j B, eqs 16 and 17 give:

Rβ SVar ) 〈φAR |φ ¯ Bβ 〉 ) cos2(θ) - sin2(θ)

Rβ ) 1, θ ) 0 and φAR and φ j Bβ are equal to ψR+ and ψβ+. When SVar Thus,

¯ β | ) |ψR ψ ¯β |φtest〉 ) |ψiRψ i + +|

R,β R,β φAR,β ) cos(θ) ψ+ + sin(θ) ψ-

(16)

R,β R,β φBR,β ) cos(θ) ψ+ - sin(θ) ψ-

(17)

R,β Unlike ηAR,β and ηBR,β where ψR,β + and ψ- have equal weights, the value of θ is determined from the variationally converged unrestricted SCF BS solution. The BS wave function is given j Bβ |, where the bar over φBβ indicates that it is as |φBS〉 ) |φAR φ occupied by a spin-down electron. The substitution of eqs j Bβ | yields 14-17 in |φAR φ

(23)

indicating a closed shell. It is not the BS wave function and is Rβ ) 0, θ ) of no further interest. On the other hand, when SVar π/4 and

¯ β | ) |ηR η¯ β | |φtest〉 ) |ψiRψ i A B

They are located on different A and B sites of the molecule respectively. The sites A and B correspond to the two BVD verdazyl rings. In the BS case, eqs 14 and 15 do not apply. However, one may define an equivalent set of similar BS orbitals as51

(22)

(24)

j βi coincide with the regular magnetic orbitals, Thus, ψRi and ψ j βB, indicating that they are not the BS magnetic orbitals. ηRA and η j βi will only correspond to φAR and φ j Bβ The orbital pairs ψRi , ψ Rβ 51 when SVar is different from 0 or 1. In this case, each unpaired electron is partially localized on a different verdazyl ring of the molecule, and the interaction between them is weakened. The BS calculations show that there is only one pair of orbitals that has an overlap that is different from 0 or 1. This pair corresponds to the unrestricted Kohn-Sham orbitals number 50 and 51. Equation 18 may be put in the equivalent form51

|ΦBS ) Cn | 1Φn + Cion | 1Φion + CT | 3Φ

(25)

with

Cn )

1 1 (cos2 θ + sin2 θ) ) √2 √2

(26)

1 [cos2 θ - sin2 θ] √2

(27)

Cion )

1 |Φ 〉 ) [cos2 θ + sin2 θ]| 1Φn〉 + √2 1 [cos2 θ - sin2 θ]| 1Φion〉 + √2sin θ cos θ| 3Φ〉 (18) √2 BS

CT ) √2sin θ cos θ

(28)

2 C2n + Cion + CT2 ) 1

(29)

and

where

1 R β [|ηAη¯ B | - |η¯ Aβ ηBR |] √2

| 1Φn〉 )

(19)

In addition, comparison of eqs 22 and 27 yields

| 1Φion〉 )

1 R β [|ηAη¯ A | + |ηBRη¯ Bβ |] √2

(20)

and

| 3Φ〉 )

1 R β [|ηAη¯ B | + |η¯ Aβ ηBR |] √2

(21)

First, one has to identify φAR and φ j Bβ in the sets of spin-up j βi }. The ψRi and and spin-down molecular orbitals, {ψRi } and {ψ β j ψi orbitals are grouped as corresponding pairs that have

Cion )

1 Rβ SVar √2

(30)

Rβ , Cn is a constant, and Thus, Cion is readily obtained from SVar CT is determined by making use of eq 29. The overlap integrals and percent character of the BS wave function, obtained from the square of the calculated coefficients as a function of φN2C3C3′N2′, are given in Table 3. Rβ 2 , Cn2, Cion The SVar , and CT2 allow the analysis of the BS wave function in terms of the conventional molecular orbitals. Table Rβ is small and tends toward 0 around φN2C3C3′N2′ 3 shows that SVar

2016



TABLE 3: BS Wavefunction Overlap, Percent Character, and Diradical Indexa φN2C3C3′N2′

Rβ SVar

C2n

2 Cion

CT2

RBS

0 10 20 30 40 50 60 70 80 90

0.157 0.143 0.107 0.061 0.009 0.042 0.087 0.123 0.145 0.153

50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00

1.30 1.00 0.60 0.20 0.04 0.09 0.38 0.75 1.05 1.25

48.77 48.98 49.42 49.81 49.99 49.91 49.62 49.25 48.95 48.84

97.537 97.962 98.847 99.626 99.992 99.827 99.239 98.491 97.892 97.671

a

2 C2n, Cion , CT2, and RBS are given as a percentage.

Figure 5. Diradical index and percent triplet character in the BVD BS wave function. Rβ When SVar ) 0, RBS ) 100%. However, when there is Rβ is unity, and RBS ) 0%, indicating no complete overlap, SVar diradical character. Table 3 also includes the RBS values over the entire range of torsional angles. From these values, the BS calculations predict that BVD has a very large diradical character ranging from 97.5 to 99.9%. This is another indicator that the spin-spin interaction between the two verdazyl rings is very small. The triplet character of the BS wave function, CT2, follows the same trend as RBS and ranges from 48.8 to 49.9%. Because 2 is very small, CT2 is almost equal to the neutral diradical Cion character Cn2, which is always 50%. A plot of CT2 and RBS is illustrated in Figure 5. As RBS increases, the triplet character, CT2, also increases at 2 . It is the expense of the antiferromagnetic ionic character, Cion 2 well-known that Cion in the BS wave function is the cause of the overestimation of the stability of the singlet state.51 Therefore, its decrease will cause JY and JN to become smaller. One is now in a position to explain why the JY, JN, and ∆EBS values are very close to one another. Because the BVD triplet state has two odd electrons, Smax 2 ) 1. Consequently, eq 2 becomes JN ) EBS - EHS ) ∆EBS. The Sˆ2 operator may be put in the more manageable form,53

Figure 3. BVD BS overlap integral and ionic character.

Sˆ2 ) Sˆ+Sˆ- - Sˆz + Sˆz2

(32)

and its expectation value becomes51,53

(

) 40°. It then increases again when φN2C3C3′N2′ increases from 2 is 40 to 90° as illustrated in Figure 3. The ionic character Cion Rβ also small and follows the same trend as SVar. At its maximum, 2 near 0 and 90°, Cion constitutes only about 1.275% of the total BS wave function. The Slater determinants of |1Φion〉 in eq 20 show that both electrons of BS magnetic orbitals are antiparallel and on the same side of the molecule. This implies a transfer of charge from one verdazyl ring to the other, creating a local ionic character as shown in Figure 4. This ionic character, albeit very small, breaks the symmetry, leads to an antiferromagnetic pathway, and favors the singlet state.51 Neese has defined the diradical index51 Rβ Rβ RBS ) 100(1 - SVar )(1 + SVar )

)(

)

NR - Nβ NR - Nβ 〈Sˆ2〉 ) + 1 + Nβ 2 2

Figure 4. Diagram of the ionic wave function.

(31)

∑ niRniβ|SiiRβ|2 i

(33) where NR and Nβ are the number of R and β electrons, respectively. For the triplet state, NR ) Nβ + 2, and 〈Sˆmax 2〉 is slightly larger than 2.0. On the other hand, for the BS state, 2 Rβ 2 〉 ) 1 - |SVar | , and eq 3 takes the form54 〈SˆBS

JY ) -

EHS - EBS 2 - (1 -

Rβ 2 |SVar |)

)

∆EBS Rβ 2 1 + |SVar |

(34)

Rβ The small SVar values in Table 3 cause the denominator to be approximately 1.0, and ∆EBS ≈ JY. The small interaction between the two unpaired electrons may be considered as a very weak form of chemical bonding which leads to low-lying excited states.51 To handle the multireference


Figure 6. Two singly occupied one-electron MOs of the planar BVD triplet ground state. Core orbitals are excluded in the orbital counting.


Figure 8. Diagram of the neutral wave function.

∆EST[CASSCF(2,2)] and ∆EST[CASSCF(6,4)] with the JY values, because whereas JY is obtained by using a BS hybrid density functional that intrinsically includes general correlation effects, the CASSCF procedure only incorporates static correlation. Nevertheless, one may interpret the relative upward shift of the ∆EST curves as follows. In a restricted Hartree-Fock procedure, the spatial wave functions, ψRi and ψβi are the same. Thus, the R and β superscripts are dropped. A two-electron closed-shell wave function takes the form

| 1Φ0〉 ) | 1Ψ1〉 )

Figure 7. RBS, JY, and CASSCF ∆EST as a function of φN2C3C3′N2′.

character induced by these low-lying excited states, one is forced to use CASSCF and MRDDCI methods. 4.3. Complete Active Space Self-Consistent Field Calculations. It is desirable to start the multireference CI calculations with CASSCF wave functions that incorporate static correlation. One starts with the smallest active space, CASSCF (2,2), that involves two electrons in two orbitals. At φN2C3C3′N2′ ) 0, they are the 2au and 2b1g singly occupied one-electron MOs, depicted in Figure 6, and form 2auX2b1g ) B1u state. Figure 7 a plots the singlet-triplet energy difference, ∆EST[CASSCF(2,2)], as φN2C3C3′N2′ is scanned from 0 to 180°. Its magnitude decreases as φN2C3C3′N2′ increases from 0 to 40° and becomes positive, indicating a ferromagnetic coupling. Subsequently, ∆EST[CASSCF(2,2)] increases in magnitude as φN2C3C3′N2′ increases from 40 to 90°. To ascertain that most of the static correlation is included, a second calculation was performed with a bigger CAS(6,4) reference space. It was constructed by adding the doubly filled 3b2g and 6b3u orbitals. ∆EST[CASSCF(6,4)] is also drawn in Figure 7 and is practically identical to ∆EST[CASSCF(2,2)]. Therefore, at this level, most of the static correlation arises from the two electrons in the two highest occupied MOs (HOMOs). This also agrees with the BS picture where |ΦBS〉 is made up of only the two unrestricted SOMO Kohn-Sham orbitals, numbers 49 (2aRu ) and β ). 50 (2b1g All four curves in Figure 7 show the same trend when φN2C3C3′N2′ increases from 0 to 180°. However, the ∆EST[CASSCF(2,2)] and ∆EST[CASSCF(6,4)] curves are shifted up compared to JY. In principle, it is inappropriate to compare

1 1 ion [| Φ 〉 + | 1Φn〉] √2

(35)

The equal weights given to ionic and neutral states in eq 35 overestimate |1Φion〉; this is schematically shown in Figure 4. However, consider for simplicity the CASSCF(2,2) case where the wave function includes an additional doubly excited j -|. It now becomes determinant, |1Ψ2〉 ) |ψ- ψ

| 1Φ0〉 ) C1 | 1Ψ1〉 + C2 | 1Ψ2〉〉

(36)

The use of eqs 14, 15, 19, and 20 gives

| 1Ψ2〉 )

1 1 ion [| Φ 〉 - | 1Φn〉] √2

and the total CASSCF wave function becomes

| 1Φ0〉 )

1 [(C1 + C2)| 1Φion〉 + (C1 - C2)| 1Φn〉] √2

(37)

During the CASSCF procedure, the C1 and C2 coefficients are adjusted variationally and have opposite signs. This additional flexibility allows the weight of |1Φn〉 to exceed that of |1Φion〉 imposed in the Hartree-Fock limit. An inspection of the schematic diagram of |1Φn〉 in Figure 8 shows that the two electrons are in different rings all the time and avoid each other. This induces static correlation and also favors |1Φn〉 over |1Φion〉. Consequently, this reduces the ionic character and causes ∆EST to decrease in magnitude and shift upward in Figure 7. The ∆EST[CASSCF(2,2)] and ∆EST[CASSCF(6,4)] upward shifts are sufficient to predict a ferromagnetic interaction when

2018



TABLE 4: MRDDCI Singlet-Triplet Energy Differences, ∆EST, as a Function of φN2C3C3′N2′ φN2C3C3′N2′ MRDDCI2(2,2) MRDDCI3(2,2) MRDDCI3(14,12)

0 10 20 30 40 50 60 70 80 90 a

-1311.9 -1192.2 -911.2 -654.5 -473.1 -442.8 -500.5 -577.5 -636.9 -662.4

-751.2 -632.2 -396.8 -206.5 -120.3 -130.5 -216.4 -354.1 -455.1 -486.1

-908.6 -667.9 -418.1 -285.0 12.1 -130.5 -156.3 -174.9 -183.5 -245.2

exp.

-887.6a

Experimental value for the corresponding TMBVD.6

φN2C3C3′N2′ ) 40°. At φN2C3C3′N2′ ) 0°, they are -387.5 and -404.6 cm-1, respectively. which are approximately half the experimental values of solid-state molecule. This underestimation of the singlet-triplet energy differences is expected because CASSCF computations do not account for dynamic correlation. However, when compared to Hartree-Fock and HDF methods, they provide good initial wave functions for CI-type calculations. 4.4. SORCI Results. The inclusion of dynamic correlation requires an MRCI method. Because BVD is relatively large and has 18 atoms, the SORCI method was used as an alternative to conventional MRCI techniques. It is specifically designed for large molecules that have more than 10 atoms.24 By starting with a small reference space generated from the two electrons in the two b1g and au SOMO orbitals, the MRDDCI2(2,2) method was used to calculate the singlet-triplet energy differences, ∆EST[MRDDCI2(2,2)]. They are listed in the second column in Table 4. The ∆EST[MRDDCI2(2,2)] values are all negative, indicating a singlet ground state for the entire range of torsion angles. At φN2C3C3′N2′ ) 0, ∆EST[MRDDCI2(2,2)] ) -1311.9 cm-1 as compared to -887.6 cm-1 for TMBVD. Thus, the dynamic correlation introduced by the MRDDCI2(2,2) method greatly overestimates the actual ∆EST. After being averaged, the natural orbitals were used as starting wave functions for a set of MRDDCI3(2,2) calculations. As mentioned previously, the MRDDCI3 method includes extra p,a 〉 and double excitations, from determinants of the type |ψi,j a,b 24 |ψi,p 〉. The resulting ∆EST[MRDDCI3(2,2)] decreased in magnitude for all the torsion angles. For example, Table 4 shows that at φN2C3C3′N2′ ) 0, ∆EST[MRDDCI2(2,2)] is -1311.9 cm-1, but ∆EST[MRDDCI3(2,2)] has decreased to -751.2 cm-1. It is now within 15.4% of the solid-state experimental value. Thus, p,a a,b 〉 and |ψi,p 〉 including AANOs and the doubly excited |ψi,j configurations in the CI expansion yields results that are closer to the experimental values. Figure 9 compares the MRDDCI2(2,2), MRDDCI3(2,2), and HDF-BS results. All three curves show the same trend. However, the ∆EST[MRDDCI2(2,2)] and ∆EST[MRDDCI3(2,2)] curves display a deeper minimum around 90° with ∆EST[MRDDCI2(2,2)] being much larger than ∆EST[MRDDCI3(2,2)], particularly at 0°. As a result of these deep minima, the JY and ∆EST[MRDDCI3(2,2)] curves cross at 63 and 118°. In summary, the ∆EST[MRDDCI2(2,2)] overestimates the solid-state experimental value by 47%, whereas ∆EST[MRDDCI3(2,2)] is closer and only 15.4% less than the experimental value. One should also note that at φN2C3C3′N2′ ) 0°, ∆EST[MRDDCI3(2,2)] and JY obtained by HDF-BS are very close to one another. To obtain further agreement with experiment, one must consider all the types of spin interaction pathways. In planar

Figure 9. MRDDCI2(2,2) and MRDDCI3(2,2) ∆EST and JY values as a function of φN2C3C3′N2′.

Figure 10. Planar BVD 3b2g and 6b3g orbitals and the corresponding 9e degenerate pair of the orthogonal form.

BVD, the OPπ system is delocalized over both verdazyl rings. The lowest valence electrons reside in σ-type orbitals followed by IPπ orbitals. The highest filled and lowest virtual (active) MOs are of the IPπ, OPπ, OPπ*, and IPπ* types. The highest virtual orbitals are IPπ* and σ*. This ordering is not necessarily maintained as one verdazyl ring is twisted relative to the other. Crossover between orbitals of different representations will occur as new types of interactions and bonding come into play. For example, the 3b2g and 6b3g orbitals are OPπ-OPπ and IPπ-IPπ respectively, as shown in Figure 10. Upon twisting the verdazyl rings relative to one another (φN2C3C3′N2′ ) 90°), they become the 9e degenerate pair which is a mixture of IPπ and OPπ character as listed in Table 1. In addition, their C3sC3′ bonds are now bonding with respect to one another. This provides additional antiferromagnetic pathways between the two verdazyl rings. Thus, the spin of the unpaired electron in one verdazyl OPπ-type orbital interacts with the spin on the electrons in the IPπ orbitals of the second verdazyl ring and vice versa. Such IPπ-OPπ and σ-OPπ interaction pathways become important particularly in the twisted D2 and D2d symmetries. The relevant configurations that represent these interaction pathways must be included in the CI wave function for a proper description of ∆EST across the entire range of φN2C3C3′N2′. Because planar BVD contains 14 bonding and antibonding π-type electrons, a minimal reference space should be constructed from these electrons. A (14,12) reference space that includes these IPπ-OPπ interactions is chosen. It is identical to that used by de Graaf and co-workers in their CASPT2 calculations.21,22 To include dynamic correlation from the σ-OPπ interaction a,b p,a pathways, the |ψi,p 〉 and |ψi,j 〉 configurations were added via the MRDDCI3 procedure. This ensures that configurations



TABLE 5: Calculated Singlet-Triplet Energy Splittings for Planar BVD method

Figure 11. Planar BVD 7b1u and 8ag orbitals and the corresponding 7b2 and 8a1 of the orthogonal form.

Figure 12. MRDDCI3 ∆EST, RBS, and JY as a function of φN2C3C3′N2′.

2 arising from excitations of the 7b1u , 7b22, 8ag2, and 8a12 orbitals, shown in Figure 11 are included in the CI wave functions. The starting wave functions for this MRDDCI3(14,12) calculation were taken from the AANOs of the previous SORCI calculations. These values are also listed in Table 4 and graphed in Figure 12. At φN2C3C3′N2′ ) 0°, ∆EST[MRDDCI3(14,12)] is in very good agreement with the ∆EST of TMBVD determined experimentally from EPR. The two values only differ by 2.3%. In addition, at φN2C3C3′N2′ ) 40°, the plot shows that the molecule is predicted to have a triplet ground state with a ferromagnetic exchange coupling of 12.0 cm-1. The main differences between ∆EST[MRDDCI3(2,2)] and ∆EST[MRDDCI3(14,12)] occur in the region of φN2C3C3′N2′ ) 40-140°. The ∆EST[MRDDCI3(14,12)] values decrease significantly when the two verdazyl rings are twisted. This is also the region where the IPπ-OPπ and σ-OPπ interaction pathways are most important, as illustrated in Figure 10. Just as ∆EST decreased in going from MRDDCI2(2,2) to MRDDCI3(2,2) as result of incorporating dynamic correlation a,b from the |Ψi,p 〉 and |Ψi,jp,a〉 excitations, the incorporation of configurations corresponding to the IPπ-OPπ and σ-OPπ interactions in the CI expansion caused ∆EST[MRDDCI3(14,12)] in the region of 40-140° to decrease even further. Figure 12 also shows that ∆EST[MRDDCI3(14,12)] and JY curves follow the same trend, with ∆EST[MRDDCI3(14,12)] being consistently smaller than JY. This complies with the theory that the HDF-BS technique overestimates the singlet state and increases the singlet-triplet energy gap. Finally, at φN2C3C3′N2′ ) 40 and 140°, ∆EST[MRDDCI3(14,12)] has the smallest value, whereas RBS is a maximum.

experimental (TMBVD) B3LYP/6-31G* MCSCF/6-31G* MCSCF/6-31G* MCQDPT//SCF/6-31G* MCQDPT//SCF/6-311G* MCQDPT//SCF/6-31G* UB3LYP(BS)a in CHCl3 UB3LYP(BS)a in Vacuo CASPT2 B1LYP (BS) MRDDCI3 a

reference space

10,10 12,12 10,10 10,10 12,12 14,12 14,12

energy (cm-1) -887.00 -664.54 -1084.3 -1084.3 -839.42 -909.37 -804.44 -721.1 -884.1 -923.36 -762.32 -908.6

energy (meV)

ref

-110 -82.39 -134.4 -134.4 -104.08 -112.75 -99.74 -89.37 -109.61 -114.49 -87.75 this work -112.66 this work

Averaged 〈J〉 values due to inter-ring rotation at 298 K.

This corroborates that at this angle, both the HDF-BS and MRDDCI3(14,12) calculations predict that the interaction between the unpaired electrons of the two radicals is a minimum. The calculations at all levels indicate that the absolute values of JY, JN, ∆EBS, and ∆EST decrease as φN2C3C3′N2′ increases from 0 to ∼40° and then increase as φN2C3C3′N2′ increases to 90°. A plausible explanation for this trend is that the antiferromagnetic pathways are proportional to the overlap between the orbitals of the two verdazyl rings. At φN2C3C3′N2′ ) 0°, the OPπ-OPπ overlap is maximum and is delocalized over the entire molecule. This overlap is proportional to cos2(φN2C3C3′N2′) and decreases as φN2C3C3′N2′ increases from 0 to 90°. Another form of overlap also exists between the OPπ and IPπ orbitals of the two rings. This type of overlap varies as cos2(90 - φN2C3C3′N2′) and is shown in the degenerate e-type orbitals of Figure 10. It increases as φN2C3C3′N2′ increases from 0 to 90°. In the region between 0 and 90°, the sum of both types of overlaps is a minimum. At this minimum (∼40°), the triplet state is stabilized, and the singlet state is destabilized relative to one another, leading to the prediction of a triplet ground state in the CASSCF(2,2), CASSCF(6,4), and MRDDCI3(14,12) cases. The BVD′s molecular symmetry dictates that an identical situation also occurs around φN2C3C3′N2′ ) 140°. A third factor affecting the magnitudes of JY, JN, ∆EBS, and ∆EST is the spin polarization caused by OPπ, IPπ, and σ interactions. Spin polarization is one of the most difficult effects to calculate in computational quantum chemistry.17 It is particularly more difficult to compute spin polarizations by using a restricted CI calculation than an unrestricted HDF technique. That is one of the reasons why a very high level of calculations, such as that of MRDDCI3(14,12), is required to get close agreement with the experimental results. Table 5 lists the previous ∆EST and JY calculations of BVD at φN2C3C3′N2′ ) 0°. It shows that Barone’s BS averaged J values 〈J〉, as a result of the ring rotations with respect to one another, are in very good agreement with the solid-state value of TMBVD.18 Our unaveraged JY, JN, and ∆EBS values are also in reasonable agreement with the TMBVD values. Thus, the BS methods provide quick and adequate results which help further guide us in ab initio calculations. The values in Table 5 indicate that multiconfiguration techniques followed by multireference methods should be employed to obtain accurate results. For example, the MCQDPT//SCF/6-311G* of Chung et al.,20 the CASPT2 of Illas et al.,21,22 and the MRDDCI3(14,12) results obtained from this work are all close to one another and to the TMBVD experimental value.

2020


One of the aims of this paper is to investigate whether a large benchmark molecule, such as BVD, can still be handled by the SORCI method. The very good agreement between the final MRDDCI3(14,12) results, obtained by the SORCI procedure, and TMBVD in the solid state indicates that size consistency is not a major problem. The difference (2.3%) between the two values at φN2C3C3′N2′ ) 0° is probably due to crystal packing effects in the solid state which are not taken into consideration in the calculations. The very good agreement between the final MRDDCI3(14,12) and CASPT2(14,12) results is also encouraging because a CASPT2 calculation that includes more than 14 electrons in the reference space becomes computationally intractable. However, the SORCI method is specifically designed to handle larger molecules24,25 and opens the way to compute the JY, JN, and ∆EST values of larger polyverdazyl compounds with MRDDCI techniques. 5. Summary and Conclusions The HDF-BS calculations predict a net antiferromagnetic interaction between the spins of the BVD two verdazyl rings, irrespective of their φN2C3C3′N2′ value. The computed JY and JN values are almost identical, and RBS is very large, ranging from 97.5 to 99.9%. This implies that these spin-spin interactions are small and may be considered as a weak form of chemical bonding that leads to low-lying excited states.51 To take into account the multireference character introduced by these lowlying excited states in a systematic fashion, one must use techniques that include both static and dynamic correlation, such as CASSCF and MRDDCI. In addition, configurations from the OPπ-OPπ, IPπ -OPπ, and σ-OPπ pathways are included in the CI expansion to ensure a balanced treatment of the spin-spin interactions. This is done by using a large reference space consisting of 14 electron in 12 orbitals and including configuraa,b p,a 〉 and |Ψi,j 〉 type in the CI wave function. tions of the |Ψi,p The variation of JY, JN, ∆EBS, and ∆EST as a function of φN2C3C3′N2′ is best interpreted as a result of OPπ-OPπ and OPπ-IPπ interactions. The OPπ-OPπ overlaps vary as cos2(φN2C3C3′N2′), whereas the OPπ-IPπ overlaps vary as cos2(90 - φN2C3C3′N2′). Their sum is a minimum around φN2C3C3′N2′ ) 40°, where the antiferromagnetic interactions are the smallest and the singlet state is least stable. Because of the BVD’s molecular symmetry, an identical situation also occurs around φN2C3C3′N2′ ) 140°. The instabilities of the singlet state in these two regions causes the CASSCF(2,2), CASSCF(6,4), and MRDDCI3(14,12) to predict a triplet state. Furthermore, in these regions, RBS is a maximum. Therefore, both the MRDDCI3(14,12) and HDF-BS results corroborate that the unpaired electron interactions around these angles are minimal. Spin polarization effects shift the entire JY, JN, ∆EBS, and ∆EST curves up and down and are also essential. They can be calculated with reasonable accuracy by unrestricted HDF and HDF-BS methods, which implicitly include exchange and correlation effects. However, to calculate spin polarization by ab initio post Hartree-Fock methods, large reference spaces and CI wave functions are needed to include the maximum amount of static and dynamic electron-electron correlation. This is another reason for performing the large MRDDCI3(14,12) calculations. The ∆EST[MRDDCI3(14,12)] at φN2C3C3′N2′ ) 0 is in excellent agreement with that of TMBVD determined experimentally from EPR and differs only by 2.3%. It indicates that size consistency is properly handled by the SORCI computations for such a large

Mattar and Dokainish benchmark molecule. In addition, ∆EST[MRDDCI3(14,12)] is consistently smaller than JY over most of the φN2C3C3′N2′ range. This confirms the theoretical assumptions that the HDF-BS technique increases the singlet-triplet energy gap and overestimates the singlet state. It is comforting to find that there is very good agreement between the final MRDDCI3(14,12) and CASPT2(14,12)21,22 results. CASPT2 calculations that include more than 14 electrons in the reference space become computationally very costly, difficult, and impractical. However, the SORCI method is specifically designed to deal with bigger molecules.24,25 Consequently, calculations of J values of larger magnetic compounds by the SORCI methodology are quite feasible. Acknowledgment. S.M. is grateful the Natural Sciences and Engineering Research Council of Canada for financial support in the form of Discovery (Operating) grants. H.D. thanks the University of New Brunswick for Graduate Student funding. We are also grateful to Professor Frank Neese for the use of his ORCA suite of programs. References and Notes (1) Lahti, M. Magnetic Properties of Organic Materials; Marcel Dekker: New York, 1999. (2) Miller, J. S.; Drillon, M. Magnetism: Molecules to Materials IV, 2003. (3) Miller, J. S., Drillon, M., Eds.; Magnetism: Molecules to Materials II: Molecule-Based Materials, 2001. (4) Miller, J. S., Drillon, M., Eds.; Magnetism: Molecules to Materials, 2001. (5) Polo, V.; Alberola, A.; Andres, J.; Anthony, J.; Pilkington, M. Phys. Chem. Chem. Phys. 2008, 10, 857–864. (6) Brook, D. J. R.; Fox, H. H.; Lynch, V.; Fox, M. A. J. Phys. Chem. 1996, 100, 2066–2071. (7) Dowd, P.; Chang, W.; Paik, Y. H. J. Am. Chem. Soc. 1986, 108, 7416–7417. (8) Dowd, P.; Chang, W.; Paik, Y. H. J. Am. Chem. Soc. 1987, 109, 5284–5285. (9) Hoshino, H.; Kimamura, K.; Iwamura, I. Chem. Phys. Lett. 1973, 20, 193. (10) Alies, F.; Luneau, D.; Laugier, J.; Rey, P. J. Phys. Chem. 1993, 97, 2922–2925. (11) Ullmann, E. F.; Boocock, D. G. B. Chem.Commun. 1969, 1161. (12) Koivisto, B. D.; Hicks, R. G. Coord. Chem. ReV. 2005, 249, 2612– 2630. (13) Chahma, M. h.; Macnamara, K.; van der Est, A.; Alberola, A.; Polo, V.; Pilkington, M. New J. Chem. 2007, 31, 1973–1978. (14) Mattar, S. M. Chem. Phys. Lett. 2005, 405, 382–388. (15) Mattar, S. M.; Sanford, J.; Goodfellow, A. D. Chem. Phys. Lett. 2006, 418, 30–35. (16) Neese, F. J. Chem. Phys. 2001, 115, 11080–11096. (17) Improta, R.; Barone, V. Chem. ReV. 2004, 104, 1231–1253. (18) Barone, V.; Bencini, A.; Ciofini, I.; Daul, C. J. Phys. Chem. A 1999, 103, 4275–4282. (19) Barone, V.; Cacelli, I.; Cimino, P.; Ferretti, A.; Monti, S.; Prampolini, G. J. Phys. Chem. ACS ASAP. (20) Chung, G.; Lee, D. Chem. Phys. Lett. 2001, 350, 339–344. (21) de Graaf, C.; Sousa, C.; de P. R. Moreira, I.; Illas, F. J. Phys. Chem. A 2001, 105, 11371–11378. (22) Moreira, I. d. P. R.; Illas, F. Phys. Chem. Chem. Phys. 2006, 8, 1645–1659. (23) Miralles, J.; Castell, O.; Caballol, R.; Malrieu, J.-P. Chem. Phys. 1993, 172, 33–43. (24) Neese, F. J. Chem. Phys. 2003, 119, 9428–9443. (25) Neese, F. Magn. Reson. Chem. 2004, 42, S187–S198. (26) Anderson, P. W. Phys. ReV. 1959, 115, 2–13. (27) Nesbet, R. K. Phys. ReV. 1960, 119, 658–662. (28) Noodleman, L. J. Chem. Phys. 1981, 74, 5737–5743. (29) Noodleman, L.; Davidson, E. R. Chem. Phys. 1986, 109, 131–143. (30) Soda, T.; Kitagawa, Y.; Onishi, T.; Takano, Y.; Shigeta, Y.; Nagao, H.; Yoshioka, Y.; Yamaguchi, K. Chem. Phys. Lett. 2000, 319, 223–230. (31) Rudra, I.; Wu, Q.; Van Voorhis, T. J. Chem. Phys. 2006, 124, 024103-024109. (32) Frisch, M.; Ragazos, I. N.; Robb, M. A.; Bernhard Schlegel, H. Chem. Phys. Lett. 1992, 189, 524–528. (33) Werner, H.-J.; Meyer, W. J. Chem. Phys. 1980, 73, 2342–2356.

6,6′-Dioxo-3,3′-biverdazyl Ground State with BS and SORCI Techniques (34) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053–5063. (35) Neese, F. ORCA - an abinitio, density functional and semiempirical program package; University of Bonn: Germany, 2007. (36) Bruna, P. J.; Peyerimhoff, S. D.; Buenker, R. J. Chem. Phys. Lett. 1980, 72, 278–284. (37) Hirsch, G.; Bruna, P. J.; Peyerimhoff, S. D.; Buenker, R. J. Chem. Phys. Lett. 1977, 52, 442–448. (38) Castell, O.; Garcı´a, V. M.; Bo, C.; Caballol, R. J. Comput. Chem. 1996, 17, 42–48. (39) Becke, A. D. J. Chem. Phys. 1997, 107, 8554. (40) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (41) Barone, V. Recent AdVances in Density Functional Methods 361; World Scientific Publishing: Singapore, 1995. (42) Neese, F.; Olbrich, G. Chem. Phys. Lett. 2002, 362, 170–178. (43) Weigend, F.; Haser, M. Theor. Chem. Acc. 1997, 97, 331–340. (44) Kelly, H. P. Phys. ReV. 1964, 134, A1450. (45) Kelly, H. P.; Sessler, A. M. Phys. ReV. 1963, 132, 2091.


(46) Neese, F.; Wennmohs, F.; Hansen, A. J. Chem. Phys. 2009, 130, 114108–114118. (47) Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1985, 115, 259– 267. (48) Ribas, V. W.; Roberto-Neto, O.; Ornellas, F. R.; Machado, F. B. C. Chem. Phys. Lett. 2008, 460, 411–416. (49) Bergman, D. L.; Laaksonen, L.; Laaksonen, A. J. Mol. Graph. Modell. 1997, 15, 301–306. (50) Laaksonen, L. J. Mol. Graphics 1992, 10, 33–34. (51) Neese, F. J. Phys. Chem. Solids 2004, 65, 781–785. (52) Mattar, S. M. Chem. Phys. Lett. 2006, 427, 438–442. (53) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to AdVanced Electronic Structure Theory, 1982. (54) Caballol, R.; Castell, O.; Illas, F.; de Moreira, I.; Malrieu, J. P. J. Phys. Chem. A 1997, 101, 7860–7866.

JP910643E

Assignment of the 6,6â²-Dioxo-3,3â²-biverdazyl Ground State by

Recommend Documents

Assignment of the 6,6â²-Dioxo-3,3â²-biverdazyl Ground State by