Chemical Physics Letters 427 (2006) 438–442 www.elsevier.com/locate/cplett
Assignment of the benzo-1,2:4,5-bis(1,3,2-dithiazolyl) ground state using the broken symmetry technique Saba M. Mattar
*
Department of Chemistry and Centre for Laser, Atomic and Molecular Sciences, University of New Brunswick, Fredericton, New Brunswick, Canada E3B6E2 Received 3 May 2006; in final form 16 June 2006 Available online 4 July 2006
Abstract The broken symmetry (BS) solutions, obtained from various electronic structure calculations, are used to assign the ground state of the benzo-1,2:4,5-bis(1,3,2-dithiazolyl) molecule (BBDTA). In every case, the exchange between the BBDTA’s two 1,3,2-dithiazolyl paramagnetic centers is found to be antiferromagnetic leading to a singlet ground state. The Hartree–Fock values of the exchange parameter, J, are the smallest while the local density functional ones are the largest. As expected, the hybrid density functionals give intermediate values. The BS wavefunctions are also analyzed and Neese’s diradical character index, RBS, is found to range from 63.5% to 99.8%. 2006 Elsevier B.V. All rights reserved.
1. Introduction There is current interest in the magnetic properties of heterocyclic radicals as molecular magnets, conducting and spintronic devices. A basic building block for these p-radicals is the 1,3,2-dithiazolyl (DTA) moiety, shown in Fig. 1, from which mono-, di- and tri-radicals etc. may be formed [1,2]. A large number of the mono-functional DTA radicals have been synthesized, studied by EPR spectroscopy and characterized by X-ray crystallography. For example, the most recent is the thiopheno-1,3,2-dithiazolyl (TDTA) [3]. In order to tailor DTA radicals to be efficient spintronic devices and molecular magnets, it is crucial to understand their aggregation mechanisms and self-assembly in the solid state [3]. However, it is just as important to comprehend how their individual electronic and magnetic properties are affected by their electronic structure and bonding. Thus, it is ideal to correlate their experimental spin Hamiltonian parameters, such as the g [4–6] and A tensor components [7], with the corresponding computed values. When a *
Fax: +1 506 453 4981. E-mail address:
[email protected].
0009-2614/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.06.093
good correlation between experiment and theory is achieved, one gets a clear idea of the radical’s electronic structure–bonding relationships, total and net spin density distributions. DTA mono-radicals have ground state wave functions that are predominantly single Slater determinants. Consequently, hybrid density functional (HDF) techniques, in conjunction with moderate basis sets, have been very successful in calculating their electronic structure, optimal geometries and A and g tensors [4–7]. The next step is to try and understand the structure and bonding properties of 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 while an antiferromagnetic one results in either a closed shell singlet or an open shell diradical. Diradicals may be categorized as two types. In the first, the two paramagnetic centers, are fixed and have a well known distance and orientation with respect to one another. In the second the two DTA radical centers, where the spins are localized, are disjointed and their relative orientations are not fixed. It is logical to start with a diradical that has a fixed spin–spin distance and a fixed orientation
S.M. Mattar / Chemical Physics Letters 427 (2006) 438–442
439
Fig. 1. The 1,3,2-dithiaza functional group enclosed in square brackets. The ‘d’ signifies that the unpaired electron is localized mainly on the nitrogen and two sulfur atoms.
such as the benzo-1,2:4,5-bis(1,3,2-dithiazolyl) (BBDTA) molecule shown in Fig. 2. The BBDTA molecule has been synthesized by various groups [8–11]. Its electron paramagnetic resonance (EPR) spectrum has been exhaustively studied and found to be relatively weak compared to the overall concentration of the sample. This suggested that the molecule mainly exists (99.0%) in an ‘EPR silent’ form. This may occur through dimerization [8,11], as shown in Fig. 3a. The observed EPR spectrum (1.0%) is made up of three components [8,11]. The first has three lines with a hyperfine splitting of 11.3 Gauss from a 14N nuclear isotope (nuclear spin I = 1). Each line is further split into another three with a smaller splitting of 0.68 Gauss. The second species is from the same radical but the hyperfine splittings were due to the 15N isotope (I = 1/2, natural abundance 0.36%) instead of the original 14N. From the line widths of the EPR spectrum in frozen glassy solutions, the magnitude of the dipolar interaction and distance between both BBDTA unpaired electrons may be estimated. Assuming point dipoles, it was concluded that the spectrum, is due to BBDTA dimers, similar to those shown in Fig. 3b
Fig. 2. The benzo-1,2:4,5-bis(1,3,2-dithiazolyl) molecule (BBDTA). In general, the two unpaired electrons are delocalized on the far left and right S2N sides and may have parallel or antiparallel spins.
Fig. 3. (a) Vertical stacking of two BBDTA molecules that produces a diamagnetic or ‘EPR silent’ species. Here, the four electrons couple to form two pairs. (b) Two BBDTA molecules that are shifted and stack to form a dimer. The two unpaired electrons that are closest together pair up. This leaves two unpaired electrons at the extreme left- and right-hand sides of the dimer.
[8,11]. The possibility of even longer oligomers was also not ruled out. In these models, the distance between the unpaired electrons at the opposite ends is long enough to cause their electron–electron exchange interaction to be negligible. Analysis of the third component revealed that it is made up of five lines separated by 0.35 Gauss [8,9,11]. Their relative intensities are 1:2:3:2:1. It was proposed that the five lines stem from ferromagnetic coupling of two radical spins at the ends of BBDT dimers or oligomers [8,11] in which the exchange is comparable with the nitrogen isotropic nuclear hyperfine coupling. The absence of any experimental evidence for a triplet DMS = ±2 transition around 1500 Gauss in the EPR spectra indicates that the monomeric BBDTA, shown in Fig. 2, has a singlet ground state [8,11]. The triplet state would have to be high enough in energy so it cannot be significantly populated at ambient temperatures. However, Barclay et al. found that the calculated BBDTA triplet and singlet electronic states were almost degenerate with an energy difference of approximately 0.5 kcal/mol [11]. The computations were performed with an effective core potential cep-31+g** basis set at the Hartree–Fock and configuration interaction (CI) levels. The computed singlet state
440
S.M. Mattar / Chemical Physics Letters 427 (2006) 438–442
was found to be slightly higher than the triplet [11]. To the best of our knowledge, the Heisenberg–Dirac–vanVleck (HDvV) spin–spin exchange parameter, J, was not calculated. The estimation of the singlet–triplet energy separation,DEST, is not a trivial task for such a large molecule and cannot be done by conventional single determinant techniques. Barclay et al. recognized this fact and accordingly used CI methods to get an estimate for this value [11]. As will be discussed in Section 3.1, accurate determination of DEST is even more difficult and requires extensive multi-reference CI calculations. This is not practical for this class of molecules if large basis sets are used. In this Letter, the broken symmetry (BS) technique, introduced by Noodleman [12,13], is used as an alternative to the CI methods to determine the ground state and J of BBDTA. In addition, the BS wavefunction is analyzed and Neese’s diradical index, RBS, is determined. 2. Computational details The BBDTA electronic structure of the triplet and broken symmetry states was calculated with the ORCA suite of programs [6] using Hartree–Fock (HF), exchange, exchange-correlation and hybrid density functionals. Calculations were done with a cluster of eight Linux computers communicating by means of a message passing interface protocol (MPICH). Since the BBDTA triplet state may be represented by a single determinant, it was appropriate to optimize the geometry of this state using HDF techniques. Once this was accomplished, the broken symmetry solution of the low spin state was then calculated using the same geometry. Barone’s EPR-II basis sets for H, C and N [14] in conjunction with the 6-31G(df,p) basis for S were used. They have been found to give very good results on DTA radicals [5,15–17]. 3. Results and discussion 3.1. Theoretical background Exchange coupling in molecules like BBDTA may be regarded as a weak form of chemical bonding between its left and right DTA moieties. This results in closely spaced levels and low-lying excited states. A system with low-lying excited states, by definition, requires treatment using elaborate multireference methods. An appropriate starting wave function must include static correlation and gives a proper description of the molecule’s electronic structure. This may be accomplished by complete active space (CAS) or restricted active space (RAS) self-consistent field (SCF) methods. Although these methods are an improvement over conventional SCF techniques, they do not account for the dynamic correlation required for accurate J calculation. Dynamic correlation may be added via multireference configuration interaction. The resulting calcula-
tions would be extremely time consuming for molecules of the size of BBDTA. Thus, technically the calculation of the exchange parameter, J, is not a trivial task. An alternative solution to the problem is to calculate J by the broken symmetry method. This technique, first devised by Noodleman, uses a variational single Slater determinant approach that does not go beyond SCF procedures [12,13]. Consequently, density and hybrid density functionals can be used to compute J. The method starts with the HDvV Hamiltonian [12] H HDvV ¼ 2JS A S B ;
ð1Þ
where SA and SB are the spins of the two electrons. They are used to define: S max ¼ S A þ S B
ð2Þ
and S min ¼ S A S B :
ð3Þ
The solutions, energies and spin expectation values for the high spin state, jSHS = Smax, MS = (SA + SB)æ, and the broken symmetry state jSBS = Smin, MS = jSA SBjæ are then determined by an unrestricted SCF procedure. In the weak coupling regime, Noodleman first estimated J as [12,13]: JN ¼
EHS EBS ; S 2max
ð4Þ
where EHS and EBS are the energy eigenvalues of the high spin and broken symmetry states, respectively. The method was later extended to the strong coupling limit by Yamaguchi et al. [18,19] and J was expressed in terms of the spin expectation values JY ¼
EHS EBS hS^2max i hS^2min i
:
ð5Þ
Eq. (5) properly reduces to Eq. (4) in the weak coupling limit. The method has been successfully used on various transition metal and nitroxide systems [20–23]. Although the BS method has been used to probe the magnetic properties of p-stacked DTA radicals in the solid state [24,25], to the best of our knowledge, it has never been used to estimate J for individual and unassociated DTA diradicals. 3.2. Numerical results The results of the BS calculations for BBDTA are presented in Table 1. All computations indicate that the total energy of the high spin triplet state is higher than the corresponding low spin broken symmetry state. Thus, the exchange coupling between the two DTA moieties of BBDTA is predicted to be antiferromagnetic. Table 1 shows that the BS Hartree–Fock (HF) calculations give the smallest JN and JY values. They are smaller than their DFT counterparts by almost one order of magnitude. This is partly due to spin contamination of the DFT and hybrid DFT BS states by excited singlet states,
S.M. Mattar / Chemical Physics Letters 427 (2006) 438–442 Table 1 Antiferromagnetic exchange coupling parametersa Method
EHS EBS
ÆS2æHS
ÆS2æBS
JN
JY
HF Xa1 Xa2 B88 BP86 PW91 PBE PBE0 B1LYP B3LYP
86.63 945.60 744.01 1233.69 878.11 903.82 933.63 457.79 433.09 496.71
2.045 2.006 2.009 2.005 2.005 2.005 2.005 2.011 2.011 2.010
1.046 0.747 0.813 0.637 0.769 0.761 0.755 0.970 0.972 0.956
86.6 945.6 744.0 1233.7 878.1 903.8 933.6 457.8 433.1 496.7
86.7 750.9 622.1 903.0 710.4 726.8 747.0 439.5 416.7 471.3
a
Energy differences and J values in cm1.
as evidenced by the deviation of their ÆS2æBS expectation values from the uncontaminated value of unity. This results in larger singlet–triplet energy differences which, in turn, lead to larger antiferromagnetic J values. The local Slater Xa exchange, for a homogenous electron gas, is the simplest form of exchange that can be included in the calculations. Theoretically a = 2/3, and when used, yields the results labeled Xa1 in Table 1. This leads to large negative JN and JY values and an ÆS2æBS = 0.747 indicating appreciable spin contamination of the BS singlet. It is well known that the theoretical a value for a homogenous electron gas is too small. For molecules, a value between 0.7 and 0.75 is more appropriate. Substituting a = 0.75 gives the Xa2 results in Table 1. The Xa2 row reveals that increasing the local exchange reduces the magnitude of the JN and JY values. Accordingly, the spin contamination decreases and ÆS2æBS increases from 0.747 to 0.813. The non-local Becke 88 exchange functional (B88) results in very large J values and the most severe ÆS2æBS spin contamination. The inclusion of the Perdew 86 correlation to the B88 exchange results in the BP86 exchange-correlation functional, which improves the situation. Table 1 shows that the J values become smaller and ÆS2æBS increases from 0.638 to 0.769. The use of the BP86, PW91 and PBE exchange-correlation functionals result in J values that are similar to one another. However, including some HF exchange via the PBE0, B1LYP and B3LYP hybrid density functionals results in much improved ÆS2æBS values (0.956–0.972) and the J values are also very close to one another. For example, the Yamaguchi J value, JY, which covers the whole exchange coupling range, differs only by 22.8 cm1 in going from the PBE0 to the B1LYP ACM hybrid functionals.
441
Neese has made use of the corresponding orbital transformation of Amos and Hall [26] such that the transformed r , are grouped in pairs of maximum similarity orbitals, w i [21]. In the case when the number of spin up and spin down electrons is the same, two situations arise. In the first, the corresponding pairs of MOs that are part of a closed shell a b have an overlap, S ab ii ¼ hwi jwi i 1: In the second, the BS valence bond like magnetic pairs will have an S ab ii that is significantly different from unity or zero. During the SCF procedure the coefficients crmi in Eq. (6) are adjusted leading to a variationally converged overlap integral [21], S var ¼ huA juB i
ð7Þ
between the pairs of BS magnetic orbitals uA and uB located on the left and right sides of the molecule, respectively. According to Neese, the BS wavefunction takes the form [21] B j ¼ cion j1 Uion i þ cn j1 Un i þ cT j3 Ui jUBS i ¼ juA u
ð8Þ
where the coefficients are obtained from the variational overlap as: pffiffiffi pffiffiffi ð9Þ cion ¼ S var = 2; cn ¼ 1= 2 and cT ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 c2ion c2n :
ð10Þ
Finally, the diradical character of the BS orbital is defined as RBS ¼ 100ð1 jS var jÞð1 þ jS var jÞ:
ð11Þ
This index is 100% if Svar = 0 and drops to 0% when Svar = 1 [21]. The cion, cn, cT and RBS parameters are important because they allow the analysis of the BS wave function in terms of conventional molecular orbitals. The j1Uionæ is schematically represented in Fig. 4. It is a superposition of two orbitals where both electrons are located on one half of the molecule. The half of the
3.3. Analysis of the broken symmetry wavefunction In the unrestricted SCF formalism the eigenvectors of the system form two sets of molecular orbitals with a and b spins. Every ith orbital may be expanded in terms of its atomic orbitals according to X wri ¼ crmi /m ; r ¼ a; b: ð6Þ m
Fig. 4. Schematic diagrams of the j1Uionæ, j1Unæ and j3Uæ components of the magnetic BS wavefunction jUBSæ.
442
S.M. Mattar / Chemical Physics Letters 427 (2006) 438–442
Table 2 Broken symmetry wavefunction components and diradical index Method
Svar
c2ion
c2n
c2T
RBS
UHF Xa1 Xa2 B88 BP86 PW91 PBE PBE0 B1LYP B3LYP
0.045 0.508 0.439 0.605 0.485 0.493 0.499 0.209 0.203 0.236
0.10 12.89 9.65 18.28 11.75 12.13 12.46 2.17 2.06 2.78
50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00 50.00
49.90 37.11 40.35 31.72 38.25 37.87 37.55 47.83 47.94 47.23
99.80 74.23 80.69 63.45 76.51 75.76 75.09 95.65 95.88 94.45
molecule that contains both electrons may be considered as negatively charged while the half which does not contain the electron is positively charged. This is the origin of its ionic nature and the ‘ion’ superscript. The column labeled c2ion in Table 2 gives the percent j1Uionæ character in the BS wave function (Eq. (8)). The HF computations result in the lowest ionic character that is less than 1.0%. On the other hand, exchange and exchange-correlation functionals have the largest ionic character, reaching 18.28% for B88. The three hybrid density functionals, PBE0, B1LYP and B3LYP have much less ionic character ranging between 2% and 3%. Increasing the contribution of the ionic configuration to the wavefunction produces an enhanced stabilization of the singlet state and hence a larger J value. The neutral singlet component, j1Unæ, has each of its two electrons on separate sides and with opposing spins (Fig. 4). In the BS formalism, it is a pure diradical, and according to Eq. (9), it constitutes 50% of the wavefunction regardless of the functional used in the calculation [21]. The triplet, j3Uæ, component is also pure diradical because the two electrons are always on opposite sides. Inspection of Fig. 4 shows that the difference between j1Unæ and j3Uæ is that j1Unæ is antisymmetric with respect to the exchange of its electrons, while j3U æ is symmetric. The j1Uionæ and j3Uæ character in the BS wave function always adds up to 50%. Thus, Table 2 shows that the HF computations result in the highest j3Uæ character (49.90%). The lowest triplet character stems from the BP88 functional. The PBE0, B1LYP and B3LYP hybrid functionals all yield a triplet character of 47–48%. The RBS derived by Neese reflects the total j1Unæ and 3 j Uæ character. It follows the same trend as the BBDTA triplet character. The B88 has the lowest RBS while HF gives the highest value. If one assumes that the PBE0, B1LYP and B3LYP give the most reasonable results then it is concluded that the BBDTA molecule has 96% diradical character with JY 450 cm1 and an antiferromagnetic ground state. This is in accord with the conclusions reached by Dorman et al. [8] and Barclay et al. [11] from their experimental results.
Finally, inspection of Tables 1 and 2 shows that when the exchange coupling is weak, such as is the case of the HF computations, the values of JN are almost identical to JY. This confirms the fact that Yamaguchi’s expression will become equivalent to Noodleman’s expression in the weak coupling limit. However, when functionals that emphasize the coupling are used then the two values deviate from one another with JN > JY. Acknowledgements I thank the Natural Sciences and Engineering Research Council of Canada for financial support in the form of discovery (operating) Grants. I am also grateful to Professor Frank Neese for the use of his ORCA suite of programs. References [1] A.W. Cordes, R.C. Haddon, R.T. Oakley, Adv. Mater. 6 (1994) 798. [2] J.M. Rawson, A.J. Banister, I. Lavender, Adv. Heterocycl. Chem. 62 (1995) 137. [3] A. Alberola, R.D. Farley, S.M. Humphrey, G.D. McManus, D.M. Murphy, J.M. Rawson, Dalton Trans. (2005) 3838. [4] S.M. Mattar, Chem. Phys. Lett. 405 (2005) 382. [5] S.M. Mattar, J. Sanford, A.D. Goodfellow, Chem. Phys. Lett. 418 (2006) 30. [6] F. Neese, J. Chem. Phys. 115 (2001) 11080. [7] R. Improta, V. Barone, Chem. Rev. 104 (2004) 1231. [8] E. Dormann, M.J. Nowak, K.A. Williams, R.O. Angus, F. Wudl, J. Am. Chem. Soc. 109 (1987) 2594. [9] K.A. Williams, M.J. Nowak, E. Dormann, F. Wudl, Synth. Met. 14 (1986) 233. [10] G. Wolmershaeuser, M. Schnauber, T. Wilhelm, L.H. Sutcliffe, Synth. Met. 14 (1986) 239. [11] T.M. Barclay et al., J. Am. Chem. Soc. 119 (1997) 2633. [12] L. Noodleman, J. Chem. Phys. 74 (1981) 5737. [13] L. Noodleman, E.R. Davidson, Chem. Phys. 109 (1986) 131. [14] V. Barone, Recent Advances in Density Functional Methods, World Scientific Publishing, Singapore, 1995. [15] S.M. Mattar, Chem. Phys. Lett. 300 (1999) 545. [16] S.M. Mattar, Chem. Phys. Lett. 405 (2005) 382. [17] S.M. Mattar, J. Sanford, A.D. Goodfellow, Chem. Phys. Lett. 418 (2006) 30. [18] T. Soda et al., Chem. Phys. Lett. 319 (2000) 223. [19] K. Yamaguchi, Y. Takahara, T. Fueno, in: V.H. Smith (Ed.), App. Quantum Chem., Reidel, Dordrecht, 1986, p. 155. [20] C. Adamo, V. Barone, A. Bencini, F. Totti, I. Ciofini, Inorg. Chem. 38 (1999) 1996. [21] F. Neese, J. Phys. Chem. Solids 65 (2004) 781. [22] S. Sinnecker, F. Neese, L. Noodleman, W. Lubitz, J. Am. Chem. Soc. 126 (2004) 2613. [23] R. Improta, K.N. Kudin, G.E. Scuseria, V. Barone, J. Am. Chem. Soc. 124 (2002) 113. [24] G.D. McManus et al., J. Mater. Chem. 11 (2001) 1992. [25] J. Luzon, J. Campo, F. Palacio, G.J. McIntyre, J.M. Rawson, Polyhedron 24 (2005) 2579. [26] A.T. Amos, G.G. Hall, Proc. Roy. Soc. Ser. A 263 (1961) 483.