Can Unrestricted Density-Functional Theory Describe Open Shell Singlet Biradicals?

Abstract

Unrestricted density functional theory (UDFT) can be used for the description of open-shell singlet (OSS) biradicals provided a number of precautions are considered. Biradicals that require a two-determinantal wave function (e.g. OSS state of carbenes) cannot be described by UDFT for principal reasons. However, if the overlap between the open-shell orbitals is small (the single electrons are located at different atomic centers) errors become small and, then, the principal failure of UDFT in these cases is not apparent and may even be disguised by the fact that UDFT has the advantage of describing spin polarization better than any restricted open shell DFT method. In the case of OSS biradicals with two- or multiconfigurational character (but a one-determinantal form of the leading configuration), reasonable results can be obtained by broken-symmetry (BS)-UDFT, however in each case this has to be checked. In no case is it reasonable to lower the symmetry of a molecule to get a suitable UDFT description. Hybrid functionals such as B3LYP perform better than pure DFT functionals in BS-UDFT calculations because the former reduce the self-interaction error of DFT exchange functionals, which mimics unspecified static electron correlation effects, so that the inclusion of specific static electron correlation effects via the form of the wavefunction becomes more effective.

1 Introduction

Biradicals such as p-didehydrobenzene [1] or α,3-didehydrotoluene [2, 3] have attracted considerable attention in the last ten years because they are the intermediates of the Bergman reaction [4] of enediynes and determine the biological activity of the latter compounds, in particular their potential as antitumor drugs. Their ground state is a singlet rather than a triplet [1, 2, 3, 6, 7], which seems to tune their reactivity in such a way that DNA is attacked by these biradicals in a regioselective fashion [8]. Since singlet open-shell biradicals are difficult to detect and to investigate by experimental means [1, 2, 3, 6, 7, 8, 9], most of their properties have been determined by quantum chemical calculations, which have become an indispensable tool in this connection [1, 3, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. Beside high-quality calculations carried out with wave function theory (WFT) using Coupled Cluster methods such as CCSD(T) [24] or BD(T) [25], work for the larger enediyne systems is increasingly carried out with density functional theory (DFT) [26, 27]. DFT has become an attractive tool for the investigation of enediyne biradicals because DFT and in particular unrestricted DFT (UDFT) lead to reasonable descriptions of biradicals [1, 11, 12, 13, 17, 18, 19, 20, 21, 22, 23] and in this way provides a convenient way of circumventing the use of much more expensive ab initio methods in these cases.

Actually, singlet open-shell biradicals can possess multideterminant (multireference) character, which implies that single determinant methods such as Hartree-Fock (HF) theory or correlation-corrected methods based on HF will fail to describe these molecules in an appropriate way. Reliable descriptions can only be obtained if a multireference approach such as MRCI or MRCC [28] is used or alternatively a single-determinant method covering a high amount of dynamic electron correlation is applied.  Methods such as CCSD(T) or BD(T) will partially compensate for the deficiencies of a single-determinant starting wave function and in this way provide a reasonable account of the properties of biradicals [10, 12, 16, 18, 29].

In view of these considerations it is astonishing that a single-determinant method such as Kohn-Sham DFT [26] has provided reasonable descriptions for some of the biradicals investigated. Cremer, Kraka, and co-workers obtained reasonable energies, geometries, and one-electron properties for m- and p-didehydrobenzene biradicals using UDFT [1, 3, 11, 12, 13, 29]. This observation is in line with the experience made by other authors when using UDFT for the description of open shell singlet (OSS) systems [18, 19, 20, 21, 22, 30, 31]. Gräfenstein and co-workers [11] analyzed the performance of UDFT in the case of p-didehydrobenzene in detail and showed that UDFT actually covers some static electron correlation in spite of its single-determinant nature. The success of DFT in describing singlet biradicals has led other authors to extend the application of this method also to systems where DFT (either unrestricted or restricted, UDFT or RDFT) should fail for principal reasons. Such attempts have have been the basis for confusing statements about the applicability of DFT. In this connection, those DFT investigations are particularly problematic that seem to lead to correct results. They are used to declare UDFT the method of choice for all chemical problems involving biradicals.

In this work, we will critically discuss the question to which extent DFT and in particular UDFT can be used for a reliable description of biradicals with multireference character. For this purpose, we will consider simple biradicals such as methylene in its OSS biradical state (¹B₁, 1), α, 3-didehydrotoluene (¹A″, 2) or 1, 4-didehydrobenzene (¹A_g, 3, Scheme 1). These electronic systems are typical representatives for the large classes of carbenes (silylenes, etc.) and didehydroannulenes, which are currently under intense investigations either to synthesize carbenes with increased life-times and special properties for synthesis or to design new antitumor drugs. The size of the biradicals experimentally interesting (e.g. diarylcarbenes or the p-didehydrobenzenes generated from naturally occurring enediynes such as dynemicin, calicheamicin, etc) exclude quantum chemical descriptions with high level ab initio methods and, therefore, DFT is amply used in these cases mostly without the chance of verifying DFT results by other more accurate methods. We have chosen molecules 1, 2, and 3 as suitable examples for the experimentally relevant biradicals because DFT results for 1, 2, and 3 can be compared with available experimental data and the results of high level ab initio calculations. In this work, we seek some basic answers on the performance of DFT in the case of biradicals with multireference character rather than presenting DFT descriptions for a large number of interesting systems.

According to a general scheme of classifying multireference systems suggested by Gräfenstein and Cremer [32] one can distinguish between type 0, type I, type II, and type III systems (Figure 1). A closed shell molecule without any multireference character such as benzene (type 0 systems) can be represented by a single determinant Φ₀ (in Figure 1, the 0, 1, etc. denote Φ₀, Φ₁, etc.) with weight factor w = 1 if just the Hartree potential for the electron interactions is used within Kohn-Sham DFT. By switching on the real electron-electron interaction with the help of the parameter λ (λ increasing from 0 to 1, adiabatic connection scheme [33], Figure 1), configuration state functions Ψ₁, Ψ₂, etc. corresponding to excited states are mixed into the wave function of the type 0 system thus describing dynamic electron correlation. For most type 0 systems it is sufficient to express excited configuration by just one determinant (Ψ₁ = Φ₁, Ψ₂ = Φ₂, etc.), in which an occupied orbital is replaced by an unoccupied orbital with additional nodal planes so that left-right, in-out, angular, and other correlation effects can be described. Since these higher configuration state functions possess weight factors of small but finite magnitude, the weight factor of Ψ₀ = Φ₀ becomes smaller than 1, but Ψ₀ still dominates the total wave function. For λ = 0 → 1, the energy decreases in the way that electron correlation is covered by an increasingly admixture of excited configurations Ψ₁ = Φ₁, Ψ₂ = Φ₂, etc. to the ground state wave function Ψ₀ = Φ₀ (Figure 1). The difference E(λ = 1) − E(λ = 0) is the correlation energy E_C, i.e. the energy change due to electron correlation effects (Figure 1).

Type II systems are multireference systems with strong static (non-dynamic) electron correlation [32]. This is reflected by the fact that for λ = λ_x an (or several) excited configuration state function(s) start(s) to make (a) large contribution(s) to the true wave function, which is reflected by a large increase in the corresponding weight factor(s) (Figure 1). The function Ψ₀ does no longer dominate the wave function. The strong mixing in of a higher configuration is reflected by a strong decrease of the energy at λ_x. The difference E_C covers now both the dynamic and the static electron correlation energy.

Ozone is a typical example for a type II system, because the description of its ground state wave functions requires several configuration state functions Ψ_i. RDFT as a single determinant method is no longer able to describe such a molecule correctly, which is reflected by the fact that the RDFT energy is considerably above the true energy as indicated in Figure 1. The UDFT energy, however, provides often a reasonable approximation for the true energy (Figure 1) because it can cover some multi-configurational effects. Since there is no guarantee that this is really the case for any type II system, it is the goal of this work to clarify for which type II systems UDFT represents a reasonable approximation.

There is a continuous transition from type 0 to type II systems, which implies that one can encounter molecules with weak, medium or strong static electron correlation effects. Type I systems do not take an intermediate position because their wave function is of different nature than that of either type 0 or type II systems [32]. The ground state configuration Ψ₀ still dominates the true wave function for λ = 0 → 1 as in the case of type 0 systems, however Ψ₀ is already for λ = 0 a two-determinantal wave function where the two determinants Φ₀ and Φ′₀ contribute with equal weight factors w = w′. Typical representatives of type I systems are 1 in its ¹B₁ OSS or 2 in its ¹A″ OSS state. Hence, if UDFT is able to describe a type II system such as 3 in a reasonable way, this does not automatically imply that a reasonable UDFT description is also obtained for type I systems.

The necessity of a two-determinantal wave function in the case of type I systems does not imply any long-range (static) electron correlation; on the contrary type I systems such as the ¹B₁ OSS state of methylene are molecules with just dynamic electron correlation, i.e. excited configurations Ψ₁ = Φ₁, Ψ₂ = Φ₂, etc. contribute to the true wave function just weakly similar as in the case of type 0 systems. Nevertheless, RDFT fails in these cases because as a single-determinantal method it cannot describe a single-configurational but two-determinantal wave function. Also, UDFT does not present a useful alternative for type I systems as is indicated in Figure 1. The fact that UDFT can describe a type II system such as 3 in a reasonable way does not automatically imply that a reasonable UDFT description is obtained for type I systems.

Similarly as in the case of type 0 systems, type I systems can adopt multiconfigurational character and one has to consider that Ψ₀ and excited configuration state functions Ψ₁, Ψ₂, etc. substantially contribute to the true wave function where each configuration state function is built up by two or more determinants with equal weight. Electron systems that require such a multi-configurational multi-determinantal wave function are denoted as type III systems (Figure 1). Biradical 2 is a representative of a type III system and as such the most difficult biradical to be described. Again, there is a continuous transition from type I to type III systems. In view of the particular nature of the type I and type III wave functions, we have to clarify whether DFT can describe type I and type III systems at all.

It is noteworthy that for the classification of atoms and molecules into type 0, I, II, ands III systems elements and techniques of WFT are used and connected with the DFT description. In this connection, however it is important to realize that electron correlation is differently defined in WFT and DFT. In the former case, it is defined as the difference between the exact energy of the interacting many-electron system and the energy of a suitable reference system, which is normally given by the HF description. Both WFT and Kohn-Sham DFT use the energy of a single-determinantal wave function as reference. This reference function is optimized with respect to the total energy for HF while in KS theory the reference function is optimized to provide the best ρ(r), which means that the reference energy is more positive for KS than for HF; hence the exact KS correlation energy is larger (by its absolute value) than the exact correlation energy in WFT. Also it has to be considered that in DFT the kinetic energy is just calculated for the non-interacting electrons. Hence the total correlation energy covers a kinetic energy correction due to the correlated movement of the electrons.

2 Ways of Describing Biradicals by Unrestricted Hartree-Fock Theory

Since the procedures used to describe biradicals by UDFT originate in HF theory [34], we will briefly review how UHF is used in this connection. The initial guess for an OSS state is normally chosen to be of closed-shell nature. The spin symmetry must be broken explicitly if the singlet biradical is to be described by UHF. This is done by manipulating orbitals ϕ_R = R and ϕ_S = S that describe the biradical character of the system. As a rule, R and S will be HOMO and LUMO of the closed-shell guess. They are shown schematically for 1 and 2 in Figure 2a.

The unrestricted wave function takes the following form:

(1)

where $\mathcal{A}$ is the antisymmetrizer, ${\Psi }_{\mathrm{core}}^{\text{U}}$ denotes the (closed-shell) core part, and ${\Psi }_{\mathrm{open}}^{\text{U}}$ describes the open-shell part. We use the symbol U rather than UHF to indicate that Eq. (1) and the following equations are equally valid for UHF and UDFT. Actually, ${\Psi }_{\mathrm{core}}^{\text{U}}$ will no longer represent an exact closed-shell function in the final unrestricted wave function due to the independent optimization of the α and β spin orbitals. There are two ways to build ${\Psi }_{\mathrm{open}}^{\text{U}}$ in the initial guess. If the HOMO and LUMO are of σ and π character or, more generally, belong to different irreducible representations of the molecular symmetry, these two ways will lead to qualitatively different results.

Permuted-Orbitals (PO) UDFT. The ordering of the orbitals is changed for one of the spin orientations; typically, the HOMO and LUMO are exchanged for β spin orientation (see Figure 2a and Figure 2b):

(2)

where non-overlined and overlined symbols denote α and β spin orbitals, respectively. In the special case of molecules 1 and 2, Eq. (2) takes the form

(3)

As shown in Figure 2b for 1, the resulting initial state is a mixture of a σπ singlet (S) and a σπ triplet (T) component with equal weights,

(4)

(5a)

(5b)

i.e. the spin symmetry of the closed-shell initial state is broken. Regarding spatial symmetry, the PO-UHF (PO-UDFT) reference state will belong to a one-dimensional irreducible representation, which is antisymmetric with respect to the mirror plane (^1,3B₁ for 1, ^1,3A″ for 2), in distinction to the closed-shell initial guess, which belongs to the totally symmetric irreducible representation (¹A₁ for 1, ¹A′ for 2).

Broken-Symmetry (BS) UDFT. In BS-UDFT, the initial guess is constructed as follows:

(6a)

(6b)

(6c)

with a rotation angle θ that is optimized during the SCF iterations. After the optimization, ϕ_r and ϕ_s are the localized counterparts of ϕ_R and ϕ_S, respectively (see 3 in Figure 2a). They resemble the Generalized-Valence-Bond (GVB) orbitals as is shown in Figure 2b. The reference state ${\Psi }_{\mathrm{open}}^{\text{BS−U}}$ can be rewritten as

(7)

The resulting state is a mixture of the |$\sigma \overline{\sigma }$〉 and |$\pi \overline{\pi }$〉 singlet states and the|$\sigma \overline{\pi }$〉^T (M_S = 0) triplet state, i.e. the spin symmetry is broken as in the case of PO-UHF (PO-UDFT). With respect to spatial symmetry, the BS-UHF (BS-UDFT) reference state does not belong to any irreducible representation but is part of a reducible representation that consists of a ¹A₁ and a ³B₁ contribution for 1 and a ¹A′ and a ³A″ contribution for 2. However, it is important to note that the BS-UHF (BS-UDFT) reference function is not completely asymmetric but represents an irreducible representation of a mixed spin-space symmetry group where all reflections at the mirror plane are combined with a simultaneous flip of all spins in the system. The PO-UHF (PO-UDFT) wave function is not symmetric with respect to this combined symmetry group.

The BS-UHF (BS-UDFT) wave function (7) covers in principle the GVB wave function and, therefore, is able to describe two-configurational problems such as biradicals where the rotational angle θ determines the amount of biradical character. Since for a given problem it is sometimes difficult to determine the optimized θ value and the exact form of (7), we use the natural orbital occupation numbers (NOON) of the mixing orbitals to calculate the biradical character covered by the BS-UHF (BS-UDFT) wave function.

Both for PO-UHF (PO-UDFT) and for BS-UHF (BS-UDFT), the wave function belongs to an irreducible representation of an appropriately chosen symmetry group, and the breaking of the spin symmetry does not imply a reduction of the molecular symmetry, i.e. the molecule remains mirror-symmetric. If the biradical does not possess any symmetry as in the OSS state of 90°-rotated alkene 4, the spatial symmetry of the PO-UHF (PO-UDFT) wave function vanishes while the spin symmetry remains. The BS-UHF (BS-UDFT) wave function obtains, beside the three components of Eq. (7), a contribution of the state |$R\overline{S}$〉^S. The symmetry of the BS-UDFT wave function with respect to the combined space-spin symmetry group vanishes. Eventually, the PO-UHF (PO-UDFT) and BS-UHF (BS-UDFT) wave functions are no longer qualitatively different.

3 Description of Open Shell Singlet Biradicals by DFT: Principal Shortcomings of UDFT

Kohn-Sham DFT relates the real electronic system (atom or molecule) to a reference system of non-interacting particles that has the same electron density distribution ρ(r) as the real electronic system The total energy E of the real system can then be represented (compare with Figure 1) as

(8a)

(8b)

where E_X is the exchange energy, E_C the correlation energy, T_s the kinetic energy, V_ext the external interaction energy, and V_H the Coulomb interaction energy. In distinction to E(λ = 0) and E_X, the correlation energy E_C cannot be calculated exactly. Thus, practical DFT calculations have to use approximate expressions for E_C that rest upon the description of the correlation hole in a homogeneous or weakly inhomogeneous electron gas A similar approximation is made for E_X. Even though E_X can be calculated exactly, the approximate expressions for E_X within the exchange-correlation energy are superior in accuracy due to a cancellation of errors between E_X and E_C.

Standard Kohn-Sham DFT calculations rely on assumption (1) and assumption (2) (see e.g [26, 27]):

• The reference state is a single-determinant state. This assumption is essential for the calculation of E_X.
• The correlation hole in the real system is described reasonably well by the model correlation hole from the homogeneous or weakly inhomogeneous electron gas. This assumption is essential for the calculation of E_C.

In singlet biradicals of type I, II, or III, one or both of these assumptions are not fulfilled. Therefore, we have to investigate whether a DFT description of these systems is still possible and in particular whether these systems can be described with UDFT.

The reference wave function (i.e. the wave function for λ = 0 in the adiabatic connection scheme [33]) of a type-I singlet biradical takes the following form:

(9)

where the open-shell part Ψ_TD is now represented by two Kohn-Sham determinants Φ₀ and Φ′₀ possessing identical weight factors (Figure 1) according to Eq. (5a):

(10a)

or in the special case of molecules 1 and 2 (Figure 2b)

(10b)

Assumption (1) of standard Kohn-Sham DFT is not fulfilled for type I systems (OSS states of 1, 2 and similar biradicals) and, therefore, neither RDFT nor UDFT is suitable for describing these systems. A proper DFT description of OSS biradicals implies that the functional E_X is modified. This does not simply concern the question of what approximate functional has to be chosen but the question of the appropriate form of the E_X functional. The exact exchange energy of Ψ_TD can be written in the form

(11a)

(11b)

The exchange energy of the triplet component |$R\overline{S}$〉^T (M_S = 0) is identical to that of the triplet components |RS〉 (M_S = 1) and |RS〉 (M_S = −1), which is readily expressed in UDFT as is E_X (|$R\overline{S}$〉) for function (2) in the spin-asymmetric state. K_RS is the exchange integral between orbitals ϕ_R and ϕ_S. Hence, Eq. (11) can be used as a starting point for constructing KS schemes for type I systems, and a number of such schemes have been suggested in the literature [35, 36, 37]. The restricted open-shell singlet (ROSS) scheme by Gräfenstein, Kraka, and Cremer [35] starts from Eq. (11a). Frank and co-workers [36] developed a scheme based on Eq. (11b). The Restricted Open-shell KS (ROKS) formalism by Filatov and Shaik [37] is a generalization of the approach of Frank and co-workers [36] and makes it possible to treat type I systems as well as other low-spin open-shell states. The three methods differ slightly in the way the correlation energy is determined. However, these differences are of minor importance for the computational results and we refer to the original publications for details [35, 36, 37].

Eq. (11b) reminds of the sum rule that was used by several authors [38] to calculate the OSS energy of type I systems. However, the formalisms of Ref.s [36] and [37] do not employ independent DFT calculations for the states |RS〉 (M_S = 1 to represent |$R\overline{S}$〉^T) and |$R\overline{S}$〉 of Eq. (11b). Instead, Eq. (11b) is used to set up one KS scheme and the system is described based on a single Kohn-Sham DFT calculation.

Table 1. Energies, geometries and NOON values for the three lowest states of methylene (1). ^a)

Method	Functional	E or ∆ E	rCH	∠HCH	NOON(a1)	NOON(b1)
3B1 State
RODFT	B3LYP	−39.16668	1.077	133.9	1	1
	BLYP	−39.14139	1.084	134.4	1	1
ROKS	BLYP	−39.14214	1.083	134.6	1	1
PO-UDFT	B3LYP	−39.16874	1.077	135.0	1.000	1.000
	BLYP	−39.14396	1.083	135.5	1.000	1.000
ROHF-DFT b	LYP	−39.14608	1.07	133.9	1	1
CAS(6,6)-DFT c	CS	−39.12067	1.102	131.0	1.000	1.000
Expt. d			1.070	134.0	1	1
1A1 State
RDFT	B3LYP e	10.5 (11.8)	1.110	101.6	2	0
	BLYP e	9.4 (10.4)	1.120	100.9	2	0
REKS	BLYP	8.8	1.119	101.1	1.980	0.020
GVB-DFT b	LYP	6.3	1.110	101.6	1.986	0.014
CAS(6,6)-DFT c	CS	6.9	1.138	100.0	1.909	0.088
Expt. d		9.1	1.111	101.9	N/A	N/A
1B1 State
ROSS-DFT	B3LYP	30.4	1.072	147.5	1	1
	BLYP	29.8	1.078	148.6	1	1
ROKS	BLYP	30.8	1.077	148.6	1	1
PO-UDFT	B3LYP	11.4	1.077	138.1	1.000	1.000
	BLYP	10.1	1.083	138.3	1.000	1.000
BS-UDFT	B3LYP	5.3	1.101	114.0	1.528	0.472
	BLYP	4.4	1.083	114.3	1.525	0.475
ROHF-DFT b	LYP	25.8	1.072	147.5	1	1
CAS(6,6)-DFT c	CS	35.2	1.131	134.1	1.000	1.000
CI f		35.4	1.077	141.5	N/A	N/A

^a) Absolute energies E in hartree, relative energies ∆E in kcal mol⁻¹, bond lengths in Å, bond angles in degrees. The natural orbital occupation numbers (NOON) are given for the two highest occupied MOs. Integer values denote NOON values fixed by the method used, floating point numbers are given when other than integer values are possible. All calculations were done with Dunning’s cc-pVTZ basis set [39] unless otherwise noted. The CAS-DFT calculations were carried out with the Colle-Salvetti (CS) functional.

^b) Geometries taken from DFT. cc-pVTZ basis at DFT geometries.

^c) From Ref. [45].

^d) Experimental geometries from Ref. [49] (¹A₁) and Ref. [50] (¹A₁: 1.075 Å, 133.9 deg; and ³B₁). E(¹A₁) − E(³B₁) from Ref. [47], converted from T₀ to T_e according to Ref. [48].

^e) Energy differences in parentheses are relative to the PO-UDFT energy of the ³B₁ state.

^f) E(¹B₁) − E(³B₁) estimated by extrapolations to large basis sets, Ref. [46a]; E(¹A₁) − E(³B₁) from SOCI+Q/ANO calculations, Ref. [46b].

summarizes the results of ROSS, ROKS, and UDFT calculations for the three lowest states (³B₁, ¹A₁, and ¹B₁) of 1. Dunning’s correlation-consistent polarized valence triple-zeta basis set cc-pVTZ [39] was used in the calculations. The GGA (generalized gradient approximation) functional BLYP [40] and the hybrid functional B3LYP [41] were employed to obtain representative values. The ROSS calculations were carried out with the program package COLOGNE 99 [42], the ROKS calculations with a local version of CADPAC 5 [43], and the standard-DFT calculations with the Gaussian 98 program package [44]. Results from CAS-DFT calculations [32, 45] are also included in Table 1 to have an appropriate reference at the DFT level that includes multi- and not just two-configurational effects in a systematic manner. For each state, the calculated NOONs for the σ and π frontier orbitals of 1 are given where an integer number indicates that the NOON value is fixed by the method used.

Since there is no experimental value available for the ¹B₁–³B₁ energy splitting of 1, Bauschlicher’s [46] configuration interaction (CI) value of 33.4 kcal/mol, which is derived from a series of calculations with varying basis set, was used as a suitable reference value. Here and in the following, positive values for singlet-triplet splittings indicate a triplet ground state, negative ones a singlet ground state. For the ¹A₁–³B₁ splitting, the experimental value of McKellar and co-workers [47] corrected for vibrational effects [48] was used. Experimental [49, 50] and calculated geometries are also compared in Table 1.

DFT provides reasonable descriptions for the two lowest states of 1, which can compete with high level WFT results. For example, the BLYP and B3LYP results for the ¹A₁–³B₁ splitting obtained with restricted DFT methods are 9.4 (REKS: 8.8 kcal/mol, Table 1), and 10.5 kcal/mol, respectively, which have to be compared with an experimental value of 9.1 kcal/mol. UDFT exaggerates the splitting (Table 1, numbers in parentheses) somewhat because the triplet state is better described than the singlet state. CAS-DFT underestimates the splitting, however the value of 6.9 kcal/mol is still better than that of CASPT2 [51] and other high level calculations [48]. The performance of DFT is a direct reflection of the fact that the biradical character of the ¹A₁ state of 1 is rather small (8.7 % at GVB-DFT; 8.8 % at CAS-DFT, see Table 1). Hence, an RDFT description of this state (0 % biradical character) does not lead to a large error provided it covers sufficient dynamic electron correlation.

The situation changes when describing the OSS state ¹B₁ of 1 and the ¹B₁–³B₁ splitting. ROSS and ROKS lead to reasonable values (30.4 and 30.8 kcal/mol, Table 1), which are in line with both the CI (33.4 kcal/mol) and the CAS-DFT result (35.2 kcal/mol, Table 1). PO-UDFT, in contrast, gives a ¹B₁–³B₁ splitting of 10.1 (BLYP) and 11.4 kcal/mol (B3LYP, Table 1), respectively, i.e. the splitting is more than 20 kcal/mol off the reference value. Also, the ¹B₁ rather than the ¹A₁ state is predicted to be the lowest singlet state of 1. For BS-UDFT, the situation deteriorates even more in the way that the ¹B₁–³B₁ splitting is predicted to be 4-5 kcal/mol (Table 1). The geometry of the ¹B₁ state is now characterized by unreasonably long CH bonds (1.09 and 1.10 Å, Table 1) and a HCH angle (114°) more than 30° smaller than predicted by ROSS or ROKS. The NOON values (1.53 and 0.47, Table 1) reveal that BS-UDFT underestimates the biradical character by 53 % while PO-UDFT actually provides the correct value of 100 %.

By inspection of Eq.s (4) and (11) as well as the diagrams for type I systems in Figure 1, the dramatic failure of UDFT in the case of 1 becomes understandable. The PO-UDFT energy of the ¹B₁ state of 1 has to be too low by at least the magnitude of the exchange integral K_RS, which is 15.6 kcal/mol in the case of 1. We conclude that even a qualitatively correct PO-UDFT-description of ¹B₁-1 is not possible. A correct description of the exchange interaction between the unpaired electrons can only be provided by ROSS and ROKS (REKS) or the more general CAS-DFT method. In the case of 1, already a two-configurational WFT description (covering the two-determinantal character of the wave function for a type I system) obtained for example with ROHF for open-shell singlet state [53c] and combined with the DFT correlation energy (ROHF-DFT [53]) already leads to reasonable excitation energies (25.8 kcal/mol, Table 1; with the 6-31G(d, p) basis 29.6 kcal/mol are obtained). At this level of theory, the ¹A₁ state is described by GVB-DFT [53], which includes nondynamic electron correlation effects thus providing for this state a better description than the ROHF-DFT description can provide for the triplet state. Consequently, the singlet(¹A₁)-triplet(³B₁) splitting is strongly underestimated at this level of theory (1.2 kcal/mol, Table 1). CAS-DFT provides a more balanced account of the excitation energies of 1 although the underestimation of the ¹A₁–³B₁ splitting seems to be typical of WFT-DFT combination methods and has been attributed to the double-counting of correlation effects [45]. The failure of BS-UDFT in the description of 1 will be discussed in detail in chapter 5.

4 Why does UDFT Give Reasonable Results for Certain Open Shell Singlet Biradicals?

Table 2. Singlet and triplet energies for α-3-didehydrotoluene (2). ^a)

				NOON values
Method	Functional	E(3A″)	∆E(1A′)	3A″, a’	3A″, a”	1A′, a’	1A′, a”
ROSS	B3LYP	-270.31971	0.8	1	1	1	1
	BLYP	−270.19931	0.9	1	1	1	1
REKS	BLYP	−270.20613	0.9	1	1	1	1
PO-UDFT	B3LYP	−270.32322	−1.0	1.000	1.000	1.000	1.000
	BLYP	−270.20152	−0.5	1.000	1.000	1.000	1.000
BS-UDFT	B3LYP	−270.32322	22.5	1.000	1.000	1.423	0.577
GVB-DFT b	LYP	−270.11665	−7.4	1	1	1	1
CAS(8,8)-DFT c	CS	−270.03909	−2.6	1	1	1	1
DDCI2 d			−2.0
Expt. e			−5... 0	1	1	1	1

^a) Absolute energies E in Hartree, relative energies ∆E in kcal mol⁻¹. All calculations were done with Dunning’s cc-pVTZ basis set [39] unless otherwise noted. The method is given for the description of the OSS state while for the triplet state either RODFT or UDFT was used. The natural orbital occupation numbers (NOON) are given for the two highest occupied MOs. Integer values denote NOON values fixed by the method used, floating point numbers are given when other than integer values are possible.

^b) Triplet calculated with ROHF-DFT. ROSS Geometries used. See Ref. [53].

^c) Ref. [32]. The Colle-Salvetti (CS) functional and the 6-31G(d,p) basis set was used.

^d) Ref. [54].

^e) Ref. [7].

lists results of DFT calculations for the lowest ¹A″ and ³A″ states of 2. The corresponding geometries are shown in Figure 3a. Experiments [3, 7] indicate that the ¹A″ state is the ground state of 2 and that the singlet-triplet splitting should be smaller than 5 kcal/mol. High-level CI calculations [21b, 54] suggest a singlet-triplet splitting of about −2 kcal/mol, i.e. the ¹A″ below the ³A″ state in line with the experimental observations.

¿From the discussion in chapter 3, one should expect that UDFT fails to describe the relative energies of the singlet and triplet states of 2. However, PO-UDFT predicts a singlet-triplet splitting of -0.5 and -1.0 kcal/mol, respectively, while both ROSS and ROKS calculations give erroneous singlet-triplet ordering as reflected by positive singlet-triplet splittings (0.8 and 0.9 kcal/mol, Table 2). Although these results seem to contradict the discussion in chapter 2, they become understandable in view of considerations (a) and (b):

(a) According to Eq. (11a), the basic failure of PO-UDFT to correctly describe the wave function of the OSS state of a type I system becomes more apparent the larger the exchange integral K_RS is. If the two unpaired electrons are well separated as in the case of 2, K_RS will become small and the basic misconception of UDFT is disguised. For 2, K_RS is equal to 0.42 kcal/mol at B3LYP/cc-pVTZ and the artificial lowering of the UDFT energy of the OSS state is negligible.

(b) It is known [32, 45, 55] that the ¹A″ and ³A″ states of 2 show spin-polarization effects in the π-electron system, which require a multi-configuration description, i.e. strictly speaking 2 is a type III rather than a type I system. In the ¹A″ state, this leads to stabilizing αα and ββ spin interactions at the single electron centers while for the ³A″ state one of these interactions is replaced by a destabilizing single-center αβ interaction. Contrary to the one-center biradical 1, biradical 2 can only be correctly described if the method used covers these spin polarization effects. ROSS and ROKS as restricted methods have to use identical α and β spin orbitals and, therefore, they are not able to cover spin polarization effects whereas UDFT can describe these effects due to the flexibility of the spin-unrestricted orbitals.

The advantages of (spin symmetry broken) PO-UDFT prove more important than the inappropriate description of E_X for the description of the singlet-triplet splitting of 2 or similar biradicals. This will no longer be true if magnetic properties or in general properties that depend on the spin density distribution are considered. A description that both accounts for static and dynamic correlations such as CAS-DFT [32, 45] is in any case preferable in such a situation because it is based on the correct wave function, describes exchange more completely, and covers also spin polarization effects. In the case of 2, CAS-DFT gives the singlet-triplet splitting of 2 as -3 kcal/mol [45] in line with difference-dedicated configuration-interaction (DDCI2) calculations [54]. GVB-DFT [53] overestimates the stability of the singlet state similar in the case of 1 (Table 2), which underlines the necessity of a balanced and, therefore, multi- rather than two-configurational approach.

5 UDFT Description of Type II Biradicals

Biradical 3 represents a type II biradical, the wave function of which can be set up with the help of Eq. (12).

(12)

The configurations |$R\overline{R}$〉 and |$S\overline{S}$〉 in Eq. (12) are non-equivalent but nearly degenerate, i.e neither |C₂| nor C₁ is small. The biradical character of wave function (12) depends on the ratio |C₂/C₁|, which becomes clear when expressing (12) in terms of the (localized) GVB orbitals ϕ_r and ϕ_s.

(13)

The first term in Eq. (13) gives the ionic terms, which become less important in the way the importance of the doubly excited state |$S\overline{S}$〉 increases. The second term describes the singlet biradical state, which gains importance with decreasing ionic terms. Hence, the ratio |C₂/C₁| is a direct indicator for the biradical character and, by this, determines the NOON values of frontier orbitals ϕ_R and ϕ_S.

The prototype of a type II system is the stretched H₂ molecule, where ϕ_R and ϕ_S are the sσ and sσ* orbitals. Type II systems possess long-range non-dynamic electron correlations, which contradicts assumption (2) of standard Kohn-Sham DFT. An RDFT treatment of these systems, although formally a correct approach, will in practice seriously underestimate electron correlation and yield too high an energy for the singlet state [11, 30] as is indicated in Figure 1. One approach to treat type II systems with DFT is to start from a two-configuration reference function that covers the static correlation effects and to use DFT for the description of dynamic correlation effects only as is done in the GVB-DFT approach [53, 56], or the two-configuration (TC)-DFT methods suggested by Borowski et al. Filatov and Shaik [58] developed the restricted-ensemble KS (REKS) method, which represents the total energy of a multi-configuration system as a weighted average over an ensemble of single-determinant states. In the case of type I systems, the REKS method becomes equivalent to the ROKS method.

BS-UDFT provides an alternative way of describing type II systems (compare, Eq.s 12 and 7). The electrons in localized orbitals ϕ_r and ϕ_s (Figure 2a) of the BS-UDFT wave function (6a) avoid each other in a similar way as in the type II reference state (12). As was discussed in detail by Gräfenstein and co-workers [11], the static correlation energy is in this way expressed as a part of the exchange energy. Hence, UDFT can lead to a reasonable description of type II biradicals provided a number of points are considered.

(a) Hybrid functional such as B3LYP in combination with BS-UDFT cover static electron correlation effects in a better way than GGA functionals such as BLYP. This becomes obvious when calculating the NOON values of the BS-UDFT description of 3: BLYP predicts a biradical character of just 51 % and by this 20 % below the corresponding B3LYP value (

Table 3. Singlet and triplet energies for 1,4-didehydrobenzene (3). ^a)

				NOON values
Method	Functional	E(3Blu)	∆E(1Ag)	3Blu, blu	3Blu, ag	1Ag, blu	1Ag, ag
RDFT	B3LYP	−230.95254	13.1	1	1	2	0
	BLYP	−230.85665	−0.1	1	1	2	0
REKS	BLYP	−230.86028	−5.7	1	1	1.723	0.277
BS-UDFT	B3LYP	−230.95254	−2.6	1.000	1.000	1.281	0.719
	BLYP	−230.85606	−4.5	1.000	1.000	1.492	0.508
BS-UDFT + SFb	B3LYP	−230.95254	−4.9
	BLYP	−230.85606	−7.3
PO-UDFT	B3LYP	−230.95254	38.9	1.000	1.000	1.000	1.000
GVB-DFT	LYP	−230.77758	−10.4	1	1	1.179	0.821
CAS(8,8)-DFTc	CS	−230.70904	−2.5	1.000	1.000	1.242	0.758
Expt. d			−3.5 ± 0.5	1	1	not known

^a) Absolute energies E in Hartree, relative energies ∆E in kcal mol⁻¹. All calculations were done with Dunning’s cc-pVTZ basis set [39] unless otherwise noted. The method is given for the description of the singlet state while for the triplet state either RODFT or UDFT was used. The natural orbital occupation numbers (NOON) are given for the two highest occupied MOs. Integer values denote NOON values fixed by the method used, floating point numbers are given when other than integer values are possible. UDFT results are from Ref. [11].

^b) BS-UDFT + sum formula.

^c) Ref. [32]. The Colle-Salvetti (CS) functional and the 6-31G(d,p) basis set was used.

^d) Ref. [6], converted to T_e according to Ref. [11].

) or the value of high level WFT calculations [10]. The low biradical character is related to the fact that BLYP leads to a more stable RDFT description of 3 than B3LYP [11]. Hence, BLYP emphasizes the closed-shell singlet structure more than B3LYP, which is also indicated by the calculated geometry of the ¹A_g state of 3 (Figure 3b): CC bond length alternation is enhanced, which is indicative of enhanced through-bond coupling of the single electrons and a reduction of the biradical character [12, 13b].

Recently, we have demonstrated with the help of electron density studies using ab initio densities as suitable reference densities that the self-interaction error of commonly used exchange functionals leads to the inclusion of nondynamic electron correlation effects in an unspecified way thus suppressing an efficient introduction of static electron correlation via a two-(multi-) configurational wave function. Since hybrid functionals such as B3LYP reduce the self-interaction error, static electron correlation effects mimicked by the exchange functional play a smaller role and, accordingly, can be more effectively introduced via the two-(multi-) configurational form of the wavefunction [59].

(b) The use of hybrid functionals leads to another advantage. Eq. (7) shows that the UDFT wave function for the OSS state is a superposition of singlet and triplet functions. For spin-unrestricted WFT, it was suggested that the influence of the triplet contribution is compensated by calculating the energy ES of the biradical state as an average using the sum formula (14) [38]

(14)

(15)

where E_T is the energy of the corresponding triplet state and ${\hat{S}}^2$ is the operator of the total spin. Application of this formula for pure DFT makes little sense, it becomes however useful as soon as hybrid functionals are used. As discussed in Ref. [11], Eq. (14) can be exploited to improve BS-UDFT results for type II systems.

(c) The magnitude of the singlet-triplet splitting provides information about the reliability of the BS-UDFT results, which are reasonable in the case of small splittings (< 3 kcal/mol), but become unreliable for large splittings (> 5 kcal/mol). Appropriate test calculations have to clarify this point in each case.

Table 3 shows RDFT, UDFT, REKS, GVB-DFT, and CAS-DFT results for the lowest singlet (¹A_g) and triplet (³B_1u) states of 3. The corresponding geometries are shown in Figure 3b. It is known from experiment [6] that 3 possesses a singlet ground state, the singlet-triplet splitting being −3.5± 0.5 kcal/mol. RDFT gives a singlet-triplet splitting of 13.1 kcal/mol for B3LYP and -0.1 kcal/mol for BLYP, respectively. Hence, the ¹A_g state is above or at equal level as the ³B_1u state contrary to experimental observations. Inspection of Figure 1 reveals that this failure of DFT results from the fact that within the single determinant approach the influence of a second (or additional) configuration(s) is not considered (NOON values of 2 and 0 for the b_1u-HOMO and the a_g-LUMO, Table 2) and, therefore, the stabilization of the ¹A_g state resulting from static electron correlation is missed. BS-UDFT partially compensates this error by adding 72 % biradical character to the DFT wave function (NOON values: 1.28 and 0.72, Table 3).

The B3LYP functional is more sensitive to indicate this basic failure of RDFT and, therefore, RDFT/B3LYP deviates more from the true description and becomes more unstable than RDFT/BLYP. This has lead to the conclusion [17, 18, 19, 20] that hybrid functionals are inappropriate for a description of 3 and related systems. Actually, if the basic error in the RDFT description is corrected by using BS-UDFT (Figure 1) it will turn out that hybrid functionals such as B3LYP combined with a sufficiently large basis set lead in general to a higher accuracy than GGA functionals such as BLYP [11]. BS-UDFT/B3LYP gives a singlet-triplet splitting of -2.6 kcal/mol, which is close to the CAS-DFT result of -2.5 kcal/mol and also close to the experimental value (Table 3) [6]. BS-UDFT/BLYP similar as REKS/BLYP exaggerate the stability of the singlet state and by this give too negative splittings (-4.5 and -5.7 kcal/mol, Table 3). REKS/BLYP itself underestimates the biradical character of the ¹A_g-state (28 %, NOON values 1.72 and 0.28, Table 3) by 30-40 % relative to high level ab initio methods [10]. Consequently, the REKS/BLYP geometry of 3 (Figure 3b) resembles more a HF than a CCSD(T) geometry [10].

Use of the sum formula can lead to even more negative singlet-triplet splittings, which suggests that application of the empirical correction has to be decided from case to case. We note in this connection that the sum formula is rather useful if BS-UDFT energies are compared with RDFT values within a reaction system [11]. The largest singlet-triplet splitting is obtained by GVB-DFT, which again is a consequence of the the inadequate description of the triplet state by ROHF-DFT.

A PO-UDFT calculation for the singlet state of 3 gives at B3LYP/cc-pVTZ an energy 38.9 kcal/mol above the triplet state, i.e. about 40 kcal/mol off the correct value. The reason is that PO-UDFT is basically inappropriate to describe type II systems: the PO-UDFT KS wave function contains no contribution of the sought type-II singlet state but describes the type I ¹B_1u (OSS) state of 3. This state has a relatively high energy because it suffers from ionic terms, which lead to large Coulomb repulsion (see Eq. 16).

(16)

Clearly, PO-UDFT must not be used for the description of type II systems. Also, the use of RDFT is erroneous in the case of a type II system with large biradical character. Of course, one could imagine a situation in which the biradical character of a given system is not known (unknown |C₂|) and therefore RDFT is applied. By checking the stability of the RDFT solution [60] it becomes obvious that UDFT has to be applied in the case of 3 and an erroneous description can be avoided in this way.

6 Use and Misuse of UDFT

The motivation for using BS-UDFT results from the necessity of describing type II systems in the most cost-efficient way. This can lead to problems and the usefulness of this approach has to be checked from case to case. By no means is it justified to derive from a successful BS-UDFT description of type II systems the validity of this description in the case of type I biradicals. Nevertheless this has been done in the literature [23] and, therefore, these attempts deserve special consideration.

Inspection of Eq. (7) reveals that ${\Psi }_{\mathrm{open}}^{\text{BS−UDFT}}$ does not contain the state of interest, i.e. |$R\overline{S}$〉^S. Instead, states are contained in the BS-UDFT wave function, for which the orbitals ϕ_R or ϕ_S are doubly occupied. It depends on the relative energies of the |$R\overline{R}$〉, |$S\overline{S}$〉, and |$R\overline{S}$〉^T states and the overlap between the |$R\overline{R}$〉 and |$S\overline{S}$〉 states (which relates to the matrix element K_RS) whether the BS-UDFT state is higher or lower in energy than the |$R\overline{S}$〉 UDFT state.

For 1, the |$R\overline{R}$〉 = |$\sigma \overline{\sigma }$〉 state is the ¹A₁ state, which is lower in energy than the ¹B₁ state to be described. This explains why BS-UDFT results in a lower energy for the ¹B₁ state, which is 5 kcal/mol below the PO-UDFT state thus increasing the UDFT error to 25 kcal/mol (Table 1). For 2, in contrast, the |$R\overline{R}$〉 and |$S\overline{S}$〉 states (Figure 2a) are higher in energy than both the |$R\overline{S}$〉^T and |$R\overline{S}$〉^S states, which is due to the higher Coulomb repulsion energy in the |$R\overline{R}$〉 and |$S\overline{S}$〉 states. As a consequence, the BS-UDFT state has a higher energy than the PO-UDFT state. One gets with the B3LYP functional a BS-UDFT singlet-triplet splitting of 22.5 kcal/mol, i.e. the energy of the BS-UDFT state is by 23.5 kcal/mol higher than that of the PO-UDFT state (Table 2). The biradical character (58 %, Table 2) is strongly underestimated, which is in line with the fact that the BS-UDFT wave function (7) gives ionic terms now too much weight.

A stability test [60] reveals that the BS-UDFT solution is internally unstable with respect to a rotation of the orbitals ϕ_r and ϕ_s (lowest eigenvalue λ of the stability matrix: -0.0636, see also Figure 4). This orbital rotation describes a breaking of the symmetry with respect to the combined spin-space symmetry group, which was discussed for the BS-UDFT state in chapter 2. This explains observations made by Schreiner and Prall [23]: If 2 and 3 are described with BS-UDFT employing an exactly planar (i.e. mirror-symmetric) starting geometry, one will obtain the BS-UDFT solution characterized by a high relative energy and a geometry that is unstable with respect to a breaking of the mirror symmetry. Upon distorting the starting geometry of 2 and 3 slightly toward non-planarity, the symmetry of the BS-UDFT initial guess is broken (Figure 4). The sought state |$R\overline{S}$〉 enters into the initial guess, and the final BS-UDFT solution takes the form of |$R\overline{S}$〉. The energy is lowered and the (correct) C_s symmetry of 2 is reduced to C₁ symmetry.

Schreiner and Prall [23] concluded that 2 may be slightly non-planar. Furthermore, they argued that pure DFT functionals are more appropriate for the description of biradicals than hybrid functionals because some of the discrepancies observed for the planar starting geometry occur only for the B3LYP functional but not for the BLYP and BPW91 functionals. The discussion of BS-UDFT shows that both of these claims are erroneous and have to be rejected. The correct choice of the UDFT initial guess (PO-UDFT rather than BS-UDFT) confirms that the geometry of 2 is planar [61]. As in the case of the RDFT description of 3, the results for B3LYP change more strongly with an incorrect use of DFT than those for pure DFT functionals, whereas DFT, when applied correctly, gives reasonable or even better results for hybrid functionals.

In summary, we state that UDFT, if correctly applied, is in many cases an appropriate method for the description of singlet biradicals while it will fail in some cases for principal reasons. It is possible to decide in advance whether a UDFT description is appropriate and, if so, which kind of UDFT approach has to be used:

• For a type I system, UDFT fails for principal reasons and should not be used unless special precautions are taken. For example, PO-UDFT can be used in some cases if the overlap between the open-shell orbitals is small, i.e. the two unpaired electrons sit at different, preferably non-neighboring atoms. In this situation the exchange integral K_RS will be negligible and the major deficiency of UDFT is disguised so that a reasonable description seems to be obtained. The PO-UDFT description of rotated alkene 4 represents such a case The use of BS-UDFT in these cases is always erroneous.
• For a type I system with considerable overlap between the open-shell orbitals, UDFT is always inappropriate. This applies e.g. to the OSS state of carbenes, where the two unpaired electrons are located both at the same atom. Correct results can only be obtained with ROSS [35] or ROKS [37].
• For type II systems, BS-UDFT will give a reasonable description provided the corresponding singlet-triplet splitting is small (< 5 kcal/mol) while PO-UDFT is inappropriate because it describes an excited OSS state of the type II system.
• Most difficult is always the description of type III systems because they possess a multiconfigurational-multideterminantal wave function. For a symmetric biradical such as 2, the PO-UDFT provides qualitatively correct energies and geometries because of a small value of K_RS (related to a small singlet-triplet splitting) and the fact that UDFT describes spin polarization better than any restricted open-shell DFT method. BS-UDFT will always be erroneous as reflected by an unreasonably high OSS energy. Also, ROSS and REKS fail in these cases because they cannot describe spin-polarization correctly.
• Neither type I, type II or type III systems can be described by RDFT. This trivial fact was overlooked for some time (see, e.g. Ref. [20]). In the case that the biradical character of a given molecule is not known one can avoid an erroneous use of RDFT by making an appropriate stability test [11, 59] of the RDFT energy.
• The calculation of NOON values helps to determine the biradical character of a molecule and to find the appropriate method for description. This procedure is regularly used when analyzing CASSCF wave functions, however its usefulness in connection with BS-UHF or BS-UDFT descriptions is hardly exploited.
• Hybrid functionals such as B3LYP perform better than pure DFT functionals when describing OSS biradicals simply because they reduce the self-interaction error of commonly used DFT exchange functionals. Hence, the inclusion of static electron correlation via the wave function is more effective than for GGA functionals [59].
• Lowering the molecular symmetry for the initial guess (e.g., breaking the mirror symmetry) may remedy an incorrect choice of the UDFT method (i.e. PO-UDFT or BS-UDFT). However, this is neither necessary nor reasonable. Firstly, both the erroneous breaking of the initial symmetry and the transition from a BS-UDFT to a PO-UDFT wave function (or vice versa) require unnecessary computational expenses. Secondly, the final geometry and electronic structure in such an optimization will still be slightly asymmetric, which may lead to qualitatively incorrect results for symmetry-dependent quantities (e.g., infrared intensities).
• The differences between PO-UDFT and BS-UDFT vanish in the case of C₁ symmetry of the biradical to be investigated. The question which of these approaches is more reasonable can only be answered by inspection of the resulting UDFT orbitals. However, it is always dangerous to consider the quality of calculated molecular properties as an appropriate criterion for the reliability of the UDFT approach chosen because this may disguise basic deficiencies of PO- or BS-UDFT.
• In any case it is advisable to use CAS-DFT as a more complete method to check the reliability of a given UDFT description. GVB-DFT is not reliable enough to fulfill this task as was demonstrated for type II and type III systems.

It is still a widespread belief that the reliability of UDFT results should be tested with the help of the expectation value 〈${\hat{S}}^2$〉. Various authors [11, 62, 63] have pointed out that 〈${\hat{S}}^2$〉 is of little use in this connection: a large value of 〈${\hat{S}}^2$〉 may be obtained in a case in which the singlet-triplet splitting is small and by this the BS-UDFT description actually rather reliable (see above). Apart from this one has to emphasize that 〈${\hat{S}}^2$〉 in DFT has no longer the same meaning as in WFT. a) The KS wave function corresponds to the situation of noninteracting electrons, which differs considerably from the many-particle state represented by the KS density. Hence, one cannot expect that 〈${\hat{S}}^2$〉 calculated from the KS Slater determinant has any physical significance. -b) The KS wave function can only describe one-particle quantities correctly while the expectation value 〈${\hat{S}}^2$〉 is a two-particle quantity. Hence, the value of 〈${\hat{S}}^2$〉 calculated with the KS Slater determinant is not correct and an expectation value for ${\hat{S}}^2$ that deviates from S(S + 1) does not necessarily indicate spin contamination.