Extent and limitations of density-functional theory in describing magnetic systems

F. Illas,1,2 I. de P. R. Moreira, 1,2 J. M. Bofill,1,3 and M. Filatov4 1Centre Especial de Recerca en Química Teòrica, Parc Científic de Barcelona, C/Baldiri Rexach 6, E-08028 Barcelona, Spain 2Departament de Química Física, Universitat de Barcelona, C/Martí i Franquès 1, E-08028 Barcelona, Spain 3Departament de Química Orgànica, Universitat de Barcelona, C/Martí i Franquès 1, E-08028 Barcelona, Spain 4Theoretical Chemistry Group, Department of Chemistry, Göteborg University, Box 460, SE-405 30 Göteborg, Sweden (Received 21 July 2004; published 29 October 2004 )

The understanding of magnetic interactions is essential to analyze and interpret, among others, neutron diffraction experiments and magnetic susceptibility measurements.Magnetic systems imply the presence of unpaired electrons and hence may appear in atoms, molecules or solids.In fact, magnetic interactions are at the hearth of molecular based magnets 1,2 or of high-T c superconductivity 3 and dominate the chemistry of radicals. 4Hence, a complete picture of electronic structure cannot avoid a detailed description of magnetic interactions.For a wide number of systems the magnetic moments are well localized on a given atom or group of atoms referred to as magnetic centers and an effective magnetic moment, S i , can be assigned to a given magnetic center i. 2,5,6 Magnetic interactions in these systems result from important electron-electron correlation effects 2,3,7 which, to provide accurate theoretical estimates of J ij values, must be appropriately described.Therefore, magnetic systems provide severe test cases for theoretical methods of electronic structure and, in particular, for density-functional theory (DFT) approaches. 8In this report the reliability of magnetic coupling constants obtained by DFT 9 is critically examined, and it is shown that the proper treatment of general open shell electronic states, where spin cannot be neglected, requires a revision of the current exchange-correlation functionals.
DFT promises an alternative approach to the exact solutions of the nonrelativistic n electron problem which are the different n electron ⌿͑r 1 s 1 ; ... ;r n s n ͒ wave functions (WFNs) satisfying the time-independent Schrödinger equation.Due to its relatively easy and efficient computational implementation, DFT is enjoying an ever-increasing popularity in Solid State Physics and Quantum Chemistry although many fundamental questions are usually bypassed.1][12][13] The energy of the n-electron system in the field of N fixed nuclei is given by where ␥ 1 ͑r 1 ͒ and ␥ 2 ͑r 1 , r 2 ͒ are the diagonal elements, ␥ 1 ͑r 1 ͒ = ␥ 1 ͑r 1 ; r 1 ͒ and ␥ 2 ͑r 1 , r 2 ͒ = ␥ 2 ͑r 1 , r 2 ; r 1 , r 2 ͒, of ␥ 1 and ␥ 2 , the usual spinless one-and two-electron density matrices, respectively. 10,12Note that the former is the one electron density ͑r͒ commonly used in electronic structure theory.It may seem that all the information necessary to calculate the energy is contained in the density matrices and one can avoid handling with ⌿͑r 1 s 1 ; ... ;r n s n ͒.However, to avoid preposterous unphysical results one still needs the knowledge of the WFN, which generates ␥ 1 and ␥ 2 . 14This can be formulated as in the following: where S n is the spin permutation group and P 21 is the permutation operator.In practice, the accurate prediction of the energy of a given system in a given electronic state still requires a reasonable estimate of ⌿͑r method, with all necessary and sufficient constraints (spin and space symmetries) to prevent the variational collapse. 14he procedure outlined above is the basis of the full configuration interaction (FCI) method.FCI provides the most accurate solution which for a Hilbert subspace of finite dimension (finite basis set) is indeed the best attainable solution and upon increasing the basis set dimension converges to the exact solution. 15,16quations (1) and (2) are also the starting point for DFT which aims to replace both ␥ 1 and ␥ 2 by ͑r͒.For the ground state, this wish is fulfilled by the Hohenberg-Kohn (HK) theorems 9 stating that the exact ground state total energy of any many-electron system is given by a universal-i.e., independent of the external potential V͑R , r͒-unknown, functional of the electron density only.The external potential and classical Coulomb electrostatic "self energy" terms in Eq. ( 2) are known functionals of ␥ 1 ͑r 1 ͒. 10 The electron kinetic energy term in Eq. (2) ͑E kin ͒ is an explicit functional of the complete ␥ 1 ͑r 1 ; r 1 Ј͒ matrix and the remaining contributions to the electron-electron term, which explicitly depend on ␥ 2 , are still unknown.In DFT, the latter and the nondiagonal part of E kin are usually included into the universal "exchangecorrelation" functional E XC ͓͔, which also depends on ͑r͒ only, and is the basis for the practical use of DFT.Notice that if E correlation ͓␥ 2 ͑r 1 , r 2 ͔͒ in Eq. ( 2) is forced to be zero, one obtains the Hartree-Fock energy (HF) whereas writing it in terms of ͑r͒ leads to the Kohn-Sham (KS) approach, provided that the nondiagonal terms of E kin and those arising from P 21 are all included in E XC ͓͔.In the HF method, the energy is obtained through a variational iterative procedure involving the nonlocal Fock operators. 13In DFT, the variational problem also has the HF mathematical structure; it can also be solved iteratively leading to the KS equations. 17The current implementation of DFT based methods differ in the particular way to model the unknown E XC ͓͔.
Since the first HK theorem states that there exists a oneto-one mapping between V͑R , r͒ and ͑r͒, it follows that ͑r͒ determines the exact nonrelativistic Hamiltonian.Hence one may, incorrectly, claim that ͑r͒ does also determine the ground state ⌿͑r 1 s 1 ; ... ;r n s n ͒.However, one must acknowledge that, even in the exact non-relativistic WFN, information regarding spin is introduced ad hoc to fulfill the Pauli principle and, in the general case, this cannot be recovered from ͑r͒.This is also acknowledged in Levy's constrained search since this procedure is manifestly forced to be only over wave functions of the correct permutation symmetry. 18owever, if one wishes to use the Levy's search to compute the energy difference between electronic states which differ essentially in the spin part one needs to ensure not only that the variational procedure leads to a ͑r͒ which derives from a certain antisymmetric ⌿ (the n-representability problem) but also that this ⌿ has a defined spin multiplicity.Otherwise, the original Levy's constrained search will only be able to find out the ground state and is useless to predict the energy difference between magnetic states.Some of the points above have been raised by several authors 10,11,19 and it has also been recognized that a proper treatment of the singlet-triplet energy difference required the introduction of symmetry dependent exchange-correlation functionals. 20Un-fortunately these well-grounded theoretical works have not attracted enough attention from the electronic structure community.In part the low impact of previous work is motivated by the lack of numerical evidence.The present work attempts to provide further arguments from both theory and careful computation points of view.
Apart from the fundamental problem of the correct introduction of spin into DFT, this work is also motivated by the erratic DFT description of the nondynamical electron correlation arising from configurational degeneracy.In particular, we are concerned with the effect of considering densities which are imposed to originate from pure spin adapted WFNs (spin restricted formalism; R-) as opposed to the standard KS approach where a spin polarized (spin unrestricted formalism; U-) single determinant is used to construct the noninteracting electron model system ͑r͒.To this end, a careful comparison between R-DFT and U-DFT calculations on magnetic problems is performed using a variety of forms for E XC ͓͔ and several WFN based approaches.For systems with two magnetic centers it has been shown that there exists a one to one mapping between spin eigenfunctions and the eigenstates of the Heisenberg Hamiltonian and also between the broken symmetry (BS) solutions obtained from U-formalisms and the eigenstates of the Ising Hamiltonian, respectively. 5,6,21,22The former are eigenfunctions of S 2 and S z whereas the latter are only eigenfuctions S z .This mapping justifies the use of the Ising model and offers a consistent scheme to extract the magnetic coupling constants through the BS approach, either in the U-HF or U-KS schemes. 5,6otice that the study of magnetic coupling in condensed matter through periodic approaches must rely in this mapping procedure.In the case of U-HF, the mapping procedure does always lead to consistent results as discussed below.However, in the case of DFT, the auxiliary nature of the KS orbitals makes this procedure at least controversial. 5,235][26] In particular, the spin-restricted ensemble-referenced KS (REKS) 26 and spin-restricted open-shell KS (ROKS) 25 formalisms furnish a general approach to treat low-spin open-shell states.These procedures show important analogies with the spin-restricted open shell HF and with the multiconfigurational selfconsistent field methods and are precisely the R-DFT methods considered in the present work.The REKS/ROKS implementation provides a description of the singlet-triplet gap which is free of the ambiguities encountered when BS solutions are used. 5,6 broad family of magnetic systems have been explored with various WFN methods as well as with different E XC ͓͔ using either the approximate BS approach or the theoretically grounded REKS formalism.The internal consistency of the REKS implementation has been checked by performing calculations in which the HF exchange is used and without the correlation functional.The results are always following the trend given by CASCI or CASSCF, as expected.Results are presented for examples representative of different families.Those are the H-He-H as model system, the 1,1Ј,5, 5Ј-tetramethyl-6,6Ј-dioxo-3,3Ј-biverdazyl radical, the ͓Cu 2 Cl 6 ͔ 2− molecular magnet and the La 2 CuO 4 high-T c superconductor parent compound. 5,27The WFN methods chosen in the present study are complete active space configuration interaction (CASCI), complete active space self consistent field (CASSCF) and difference dedicated configuration interaction (DDCI) extensively used to describe magnetic coupling in a broad family of compounds. 7,27In all cases, the minimum CAS including the two unpaired electrons in two (magnetic) orbitals has been considered; the same active space is employed in the REKS calculations.
The DFT methods include the local density approximation (LDA), and several methods based either on the generalized gradient approach (PW91-PW91, B-PW91 and B-LYP) or on HF/DFT hybrid approaches (B3-LYP). 28For further details see also Refs. 5 and 8. Tables I and II report the energy difference between the electronic states related to the magnetic order.The same standard Gaussian type orbital atomic basis sets are used for all methods. 29For the R-formalisms, this energy difference corresponds to the singlet-triplet gap and, hence, it is directly related to the Heisenberg J parameter.In the U-formalisms, the energy differences involve the ferromagnetic state (͉FM͘) and the space and spin symmetry broken solution (͉BS͘) which is generated via spin inversion on one center.In the particular case where Fock exchange is used and correlation neglected in E XC ͓͔, it has been shown 5,21 that J =2͓E͉͑BS͒͘ − E͉͑FM͔͒͘ since ͉BS͘ is a 50% mixture of singlet and triplet states.However, one can claim that DFT is an effective one electron theory. 30Hence, S 2 and spin density (in the absence of magnetic fields) is not rigorously defined in DFT and should not be treated in the same way as in WFN theory. 30Thus, one accepts the DFT energy of the ͉BS͘ solution and disregards the qualitatively incorrect S 2 and spin-density. 30In other words, the ͉BS͘ energy is regarded as an approximation to the exact singlet energy and one uses J = E͉͑BS͒͘ − E͉͑FM͒͘.In what follows, it is demonstrated that the last supposition, although perhaps numerically acceptable, 30,31 leads to contradictions with the rigorous firstprinciples point of view discussed above.Furthermore, this approximation owes its coincidental numeric accuracy to the incomplete cancellation of the electron self-repulsion (contained in the classical Coulomb term of the KS Hamiltonian) by an approximate exchange functional E X ͓͔. 31 First of all we comment on the results obtained using the highly accurate DDCI method.In all cases, the calculated J values are comparable to experimental available data and for H-He-H nicely reproduce the exact FCI values.The CASSCF and CASCI methods tend to systematically underestimate the superexchange contribution to J, because the dynamic electron-electron correlation effects are neglected. 7n some cases, as in ͓Cu 2 Cl 6 ͔ 2 , this leads to qualitatively incorrect results.The U-HF values (not shown) are close to CASCI or CASSCF results, provided that appropriate spin projection is carried out. 5This is because the symmetry breaking introduces a certain amount of nondynamical electron-correlation, 32 which is largely responsible for the superexchange contribution to the magnetic coupling parameter. 5,7Here, it is important to remark that the U-HF TABLE I. Singlet-triplet energy differences (J in cm −1 ) in linear H-He-H for various H-He distances.In the case of U methods the BS-FM value has been considered.WFN  A totally different situation appears when examining the family of DFT results.The difference between REKS and U-KS results appears to be rather erratic, not following any definite trend.This is completely contrary to the situation observed in WFN calculations.In this case, the R-HF calculations (not shown) neglect all correlation and result in a huge overestimation of direct exchange in J.However, U-HF calculations provide a fair estimate of the CASSCF results as indicated above (provided spin projection is carried out).In the case of DFT calculations, the numerical results would suggest that in some cases one needs to carry out a pertinent spin projection whereas in other cases this seems unnecessary because of the fortuitous coincidence of R-KS and U-KS results.To make it worse, this erratic behavior is not only dependent on the particular choice of E XC ͓͔ but also is system dependent.The implication of these results are very serious because they indicate that, in the absence of experimental references, DFT predictions of magnetic Heisenberg parameters obtained with the available approximate density functionals for E XC ͓͔ cannot be trusted.At first sight, one can argue that the deficiencies of DFT methods are related to the parametrization of the particular approximate functional chosen to represent E XC ͓͔.However, the strange behavior does also hold for LDA and PW91-PW91 which, are derived from first principles for the uniform or weakly inhomogeneous electron gas although the PW91 implementation of GGA is only approximate.Results in Table II show that large differences between R-KS and U-KS values are found for the Cu binuclear complex (as one would expect by comparison to the HF theory).However, these differences are largely reduced for the biverdazyl organic radical and almost disappear for the La 2 CuO 4 cluster model.
The reported results provide compelling evidence of inconsistencies of the current forms of the E XC ͓͔ and the need for a continued search for more accurate expressions.Also, it provides a clear numerical illustration of the problems of current exchange-correlation potentials, already pointed out by other authors. 10,11,19,20In this sense, the REKS/ROKS implementation of DFT avoids the formal problems of the unrestricted approaches and provides a proper representation of the electronic structure of atoms, molecules and solids in the nonrelativistic domain as in WFN methods.In conclusion, in order to use DFT to solve the exact, nonrelativistic, many-electron problem, one should impose the space and spin symmetry constraints as it is currently done in ab initio WFN theory.Bypassing this fundamental theoretical requirement would lead to inconsistent results.In the case of magnetic coupling, DFT can provide seemingly accurate numerical results but for the wrong physical reasons.
Financial support from the Spanish Ministerio de Ciencia y Tecnologia-projects BQU2002-04029-C02-01, BQU2002-00293 and the Ramon y Cajal program (I. de P. R. M.)-and in part from the Generalitat de Catalunya (projects 2001SGR-00043 and 2001SGR-00048) is fully acknowledged.
1 s 1 ; ... ;r n s n ͒.It is customary to expand ⌿͑r 1 s 1 ; ... ;r n s n ͒ in a known basis set and to find the expansion coefficients using the variational

TABLE II .
and DFT methods are ordered according to their accuracy; for each family the most accurate results first.Singlet-triplet energy differences (in cm −1 ) for the biverdazyl radical, the ͓Cu 2 Cl 6 ͔ 2− and the two center cluster model representation of La 2 CuO 4 .In the case of U-methods the BS-FM value has been considered.WFN and DFT methods are ordered as in TableI.Experimental results are discussed in Refs. 5 and 27.Biverdazyl ͓Cu 2 Cl 6 ͔ 2− La 2 CuO 4