Magnetic Exchange Couplings from Noncollinear Spin Density Functional Perturbation Theory
aa r X i v : . [ c ond - m a t . m t r l - s c i ] O c t Magnetic Exchange Couplings from Noncollinear Spin Density FunctionalPerturbation Theory
Juan E. Peralta and Veronica Barone
Department of Physics, Central Michigan University, Mt. Pleasant, MI 48859 (Dated: November 10, 2018)We propose a method for the evaluation of magnetic exchange couplings based on noncollinearspin-density functional calculations. The method employs the second derivative of the total Kohn-Sham energy of a single reference state, in contrast to approximations based on Kohn-Sham totalenergy differences. The advantage of our approach is twofold: It provides a physically motivatedpicture of the transition from a low-spin to a high-spin state, and it utilizes a perturbation scheme forthe evaluation of magnetic exchange couplings. The latter simplifies the way these parameters arepredicted using first-principles: It avoids the non-trivial search for different spin-states that needsto be carried out in energy difference methods and it opens the possibility of “black-boxifying”the extraction of exchange couplings from density functional theory calculations. We present proofof concept calculations of magnetic exchange couplings in the H–He–H model system and in anoxovanadium bimetallic complex where the results can be intuitively rationalized.
I. INTRODUCTION
Empirical models based on the Heisenberg spin Hamil-tonian are routinely utilized to describe the behaviorof a variety of magnetic systems. In most cases, thesesimple models are found to fit the experimental datavery well, provided that the parameters in the modelHamiltonian are chosen properly. The set of param-eters can include both, external parameters (tempera-ture, applied magnetic field, etc.), and internal parame-ters (magnetic exchange couplings, magnetic anisotropy,etc.). Internal parameters for a particular system can beobtained either by fitting experimental data or from first-principles electronic structure calculations by mappingtotal electronic energies to the energies of the Heisenbergspin Hamiltonian.
In particular, magnetic exchangecouplings, J , can be obtained considering the isotropicHeisenberg Hamiltonianˆ H = − X J ij ˆ S i · ˆ S j , (1)where ˆ S i and ˆ S j are the (localized) spin operators asso-ciated to each magnetic center.Perhaps the one of the most interesting manifestationof magnetism at the molecular scale can be found in com-plexes containing transition metal atoms. Many applica-tions have been suggested exploiting these molecular-sizemagnets, such as quantum computation units and high-density data storage. Due to the relatively large sizeof most complexes of interest, density functional theory(DFT) offers the most efficient alternative for model-ing the electronic structure of these systems from first-principles.
Several approaches had been proposed to extract J couplings from DFT energies. According to the spin-projected (SP) approach, the energies of a two-center complex A and B can be related to the J couplingas E LS − E HS = 4 S A S B J AB , (2) while in the non-projected (NP) approach , the energiesof a two-center complex S A and S B can be related to the J coupling as E LS − E HS = (4 S A S B + 2 S B ) J AB , (3)where S B ≤ S A . In Eqs. (2) and (3), E HS is the en-ergy of the high-spin state and E LS is the energy ofthe low-spin (broken-symmetry) state. Eqs. (2) and (3)can be straightforwardly generalized to a set of equa-tions for complexes with multiple magnetic centers. While the SP and NP methods are fairly popular, othermethods have been proposed in the literature such asNishino’s approach , the constrained-DFT approach ofRudra et al., , the Slater’s transition state method ofDai and Whangbo. and the local spin method of Clarkand Davidson. All these approaches rely on the eval-uation of the energy difference between two (for the sim-plest case of a bimetallic complex) or more states. Theevaluation of this energy difference is commonly done bycarrying out several self-consistent field calculations, onefor each different magnetic configuration. However, inmany cases, converging to the right target state could becumbersome, specially for systems containing multiplecenters with many magnetic configurations. Therefore,developing an approach that can be used in a “black-box” manner is crucial to systematically explore a largeset of complexes or complexes containing many magneticcenters.In this work, we present an approach for the evaluationof magnetic exchange couplings based on noncollinearspin density functional calculations that allows, in anal-ogy to response properties, to express the magnetic ex-change couplings as a derivative of the total electronicenergy of one single state with respect to an external pa-rameter, opening the possibility of “black-boxifying” theextraction of magnetic exchange couplings from densityfunctional theory calculations.
II. THEORY AND IMPLEMENTATIONA. Exchange Couplings as Energy Derivatives
Let us consider the effective interaction energy betweentwo magnetic centers A and B given by the isotropic clas-sical Heisenberg model (obtained by taking the expecta-tion value of Eq. (1)), E AB = − J AB S A · S B = − J AB S A S B cos θ , (4)where S A and S B are the (perfectly localized) magneticmoment vectors, J AB is the exchange coupling constant,and θ is the angle between S A and S B . From Eq. (4),one can trivially obtain J AB from the second derivativeof E AB with respect to θ at the equilibrium points, J AB = 12 S A S B (cid:18) d E AB dθ (cid:19) θ =0 = − S A S B (cid:18) d E AB dθ (cid:19) θ =180 ◦ . (5)These simple relations provide a direct path to the eval-uation of magnetic exchange couplings J AB from densityfunctional calculations if E AB in Eq. (5) is replaced bythe total Kohn-Sham (KS) energy of the system, E KS .Therefore, assuming that the electronic system dependson θ as an ideal Heisenberg model (the validity of thisassumption will be discussed in the next Section), onecan express the exchange coupling constant J AB in termsof an energy derivative as J AB = 12 S A S B (cid:18) d E KS dθ (cid:19) θ =0 (6)or J AB = − S A S B (cid:18) d E KS dθ (cid:19) θ =180 ◦ , (7)where the angle θ in the DFT framework is defined as theangle between the local magnetization vectors S A and S B . Another related method based on the Green’s func-tion formalism for crystals has been proposed by Liecht-enstein et al. B. Constraint Noncollinear Spin-DFT Calculations
To evaluate the dependence of E KS on θ , we first in-troduce two-component spinors as Kohn-Sham orbitals,Ψ i ( r ) = (cid:18) ψ ↑ i ( r ) ψ ↓ i ( r ) (cid:19) , (8)where ψ ↑ i ( r ) and ψ ↓ i ( r ) are spatial orbitals expanded in alinear combination of atomic orbitals, ψ ωi ( r ) = X µ c ωµi φ µ ( r ) ( ω = ↑ , ↓ ) . (9) The two-component spinors introduce the freedom in thespin-dependence of the KS system that allows for localrotations of the spin density characterized by θ = 0 and θ = 180 ◦ , i.e. noncollinear spin densities. Thelocal magnetization vectors S A and S B can be written as S A,B = Z d r W A,B ( r ) s ( r ) , (10)where s ( r ) = X i ∈ occ Ψ † i ( r ) σ Ψ i ( r ) (11)is the spin-density vector and W A,B ( r ) is a scalar weightfunction that determines each local magnetic site. It isimportant to recall that the magnetic centers A and B represent a group of one or more atoms. In Eq. (11), σ = ( σ x , σ y , σ z ) is the 2 × direction of the local magneti-zations S A and S B by means of the Lagrange multipli-ers technique. To this end, we construct the Lagrangianfunctional ΛΛ[ { Ψ( r ) } , λ A , λ B ] = E KS [ { Ψ( r ) } ] − λ A · ( S A × ˆ z ) − λ B · ( S B × ˆ e θ ) , (12)where ˆ e θ = sin θ ˆ x + cos θ ˆ z is a unity vector to which S B is constraint to be parallel to, and for simplicity S A has been chosen to be constraint to the z direction, asschematized in Fig. 1. Here 0 ≤ θ ≤ ◦ is consid-ered as an external parameter for which ( dE/dθ ) θ =0 =( dE/dθ ) θ =180 ◦ = 0. Eq. (12) can be readily generalizedfor the case of many magnetic centers and arbitrary unityvector directions. For the case of two magnetic centers(and for the purpose of this work), Eq. (12) does not im-ply any loss of generality. In Eq. (12), E KS [ { Ψ( r ) } ] rep-resents the KS energy of the system which is, in practice,a functional of the set of occupied KS orbitals { Ψ( r ) } .Stationary solutions of Λ for a given (fixed) θ imply d Λ d λ A = S A × ˆ z = , (13) d Λ d λ B = S B × ˆ e θ = , (14)and δ Λ δ Ψ † i ( r ) = δE KS δ Ψ † i ( r ) − h W A ( r ) λ A · ( σ × ˆ z ) + W B ( r ) λ B · ( σ × ˆ e θ ) i Ψ i ( r ) = 0 ( i ∈ occ) . (15)While Eqs. (13) and (14) restore the constraint condi-tions, Eq. (15) combined with the orthonormality condi-tion for the spinors yields a modified set of KS equations(in terms of two-component spinors) that include the twoadditional terms inside the square brackets on the left-hand side of Eq. (15), h T + V N + J + V xc − W A ( r ) λ A · ( σ × ˆ z ) −W B ( r ) λ B · ( σ × ˆ e θ ) i Ψ i ( r ) = ǫ i Ψ i ( r ) , (16)where T = − / ∇ is the kinetic energy, V N is theelectron-nuclei potential, J is the Coulomb (or Hartree)potential, and V xc is the exchange-correlation (XC) po-tential. The sum of the first four terms inside the squarebrackets is the standard KS Hamiltonian, while the twoadditional terms can be interpreted as a potential orig-inated in a torque exerted on the local magnetic mo-ments S A and S B . It should be noted that other ap-proaches had been proposed in the literature to con-straint the direction of the local mangetization in spinDFT calculations. It is important to note that since the constraint condi-tions are linear in the spin density vectors, the additionalterms in Eq. (16) depend implicitly on the orbitals onlythrough λ A and λ B , which simplifies the implementa-tion. (cid:84) S A S B z (cid:214) (cid:84) e (cid:214) FIG. 1: Schematic representation of the constraint vectorsemployed for the local rotations of the spin density.
Neglecting spin-orbit interaction, T , V N , and J arediagonal in the 2 × ψ ↑ i and ψ ↓ i is V xc . In aprevious work, we have generalized V xc for noncollinearmagnetizations, assuming that the XC energy dependson the local variables in the same manner as in thestandard collinear (spin-unrestricted) case, and impos-ing the condition for the XC energy to be invariant un-der rigid rotations of the spin density. In that work, wehave derived V xc for general energy functionals contain-ing a variety of ingredients beyond the local-spin den-sity approximation (LSDA) and the generalized-gradientapproximation (GGA), such as meta-GGAs and hy-brid density functionals. This same generalization isadopted throughout this work. Other implementationsbased on plane-waves and Gaussian-type or-bitals can be found in the literature for a variety ofapplications. The constraint vectors can be chosen without loss ofgenerality to lay in the x − z plane. Hence, as spin-orbit interaction is not included in the Hamiltonian, the two-component spinors are purely real. As a consequence, inthis scheme the orbital magnetization of the solutions isalways zero.The value of Λ at the stationary solutions for a fixed θ given by Eqs. (13)–(15) can be directly associated tothe energy of the KS system in the presence of the con-straints, E KS ( θ ). The only arbitrariness in the formula-tion is the choice of the weight function W A,B ( r ) and thefragments A and B . Since in practice it is necessary toevaluate the matrix elements of W A,B ( r ) in the atomicorbitals basis set φ ξ ( r ), it is convenient to define S A and S B using population analysis, S A,B = X µ,ν ( W A,B ) µν P µν , (17)where P µν is the spin-density matrix vector whose Carte-sian components are P xµν = P ↑↓ µν + P ↓↑ µν , (18) P yµν = i ( P ↑↓ µν − P ↓↑ µν ) , (19)and P zµν = P ↑↑ µν + P ↓↓ µν , (20)written in terms of the 2 × P ωω ′ µν = X i ∈ occ c ωµi c ω ′ ∗ νi ( ω, ω ′ = ↑ , ↓ ) . (21)For this work, we employ L¨owdin population analysis. From the expression for the atomic spin-population givenby Eq. (17), ( W A,B ) µν can be obtained as( W A,B ) µν = dS ηA,B dP ηµν ( η = x, y or z ) . (22)Thus, using the L¨owdin partitioning, it is straightfor-ward to show that the matrix elements of W A,B ( r ) canbe calculated from the atomic orbitals overlap matrix( S ) λσ = R d rφ λ ( r ) φ σ ( r ) as( W A,B ) µν = X λ ∈ A,B ( S / ) λµ ( S / ) νλ . (23)It is worth commenting that atomic spin-densities areless sensitive to the choice of the population method thanatomic densities. Since only the direction of the atomicspin-density is relevant in our method for the calcula-tion of magnetic exchange couplings, we expect even lesssensitivity with the coice of the population method.The set of spinors that simultaneously satisfy Eq. (16)and Eqs. (13) and (14) needs to be determined self-consistently since J and V xc depend on the spinors. Toobtain the KS Hamiltonian, we add the additional con-straint terms of Eq. (16) to the standard KS Hamilto-nian, initially using a guess for the Lagrange multipliers λ A and λ B . Then we determine the optimal λ A and λ B such that Eqs. (13) and (14) are satisfied, using the den-sity matrix obtained from diagonalizing the KS Hamil-tonian of the constraint system to evaluate S A and S B .This is carried out using the steepest descent method tofind the minimum of Λ as a function of λ A and λ B . Oncethe optimal values of the Lagrange multipliers are deter-mined, and provided that self-consistency is not achieved,we proceed to the next self-consistent iteration using thedensity matrix from the previous iteration. The processstops once the criteria for changes in the density matrixand total energy are met. Several consistency checks wereperformed to verify the robustness of our code. We haveimplemented this scheme in the Gaussian DevelopmentVersion program. Using this methodology, E KS ( θ ) was calculated forsmall values of θ around θ = 0 and θ = 180 ◦ and magneticexchange couplings were obtained from the quadratic co-efficient of a polynomial fit. It is important at this pointto mention that there are cases where existing approxi-mate density functional methods have difficulties in rep-resenting the LS state ( θ = 180 ◦ ). The Kohn-Shamdeterminant in these cases correspond to a “broken-symmetry” solution that mixes two or more eigenfunc-tions of the S operator. However, for the HS state( θ = 0), it is customary accepted that approximate den-sity functionals provide a reliable representation. Thus,even though for comparison purposes in the next Sectionwe show our results using both the HS and LS statesas reference, for the practical extraction of magnetic ex-change couplings in the DFT framework, this method isexpected to work more reliably using only the HS as thereference state. III. PROOF OF CONCEPT CALCULATIONS
We first tested our methodology in the H–He–H linearmodel system with a distance H–He of 1.625˚A, consid-ering the outer H atoms as magnetic centers A and B ( S A = S B = 1 / d E KS /dθ using the θ = 0 ( h S z i = 1) state and the θ = 180 ◦ ( h S z i = 0) state, J HS and J LS , respectively. All cal-culations were carried out with the 6-311G** Gaussianbasis set. For comparison, in Table I we include resultsfor the LSDA (Dirac exchange plus the parametrizationof Wosko, Wilk, and Nusair for correlation), the BLYPrealization of the GGA (Becke’s 1988 functional for ex-change and the correlation functional of Lee, Yang, andParr ), and for the B3LYP hybrid functional.For all functionals, exchange couplings calculated fromthe energy derivatives, J HS and J LS , are in very closeagreement to the exchange coupling calculated from theenergy difference, J ∆ E . The difference can be attributedto both, the intrinsic accuracy of the numerical differ-entiation method and to the fact that J HS and J LS areexpected to be identical to J ∆ E only in the case where DFT describes the electronic system as an ideal Heisen-berg model. The small discrepancy between J HS , J LS and J ∆ E can be understood in terms of the localized na-ture of the magnetization on the H atoms in this modelsystem and provides a measure of how well the electronicsystem mimics the behavior of an ideal Heisenberg model.In Table I we report J ∆ E as calculated from the SPformula, Eq. (2), since it offers a direct comparison be-tween J HS and J LS , and J ∆ E . It is worth mentioningthat in the ideal case of a perfect Heisenberg system,∆ E = E ( θ = 180 ◦ ) − E ( θ = 0) is related to J HS and J LS according to∆ E = 2 (cid:18) d Edθ (cid:19) θ =0 = − (cid:18) d Edθ (cid:19) θ =180 ◦ . (24)Therefore, d E/dθ provides a measure of ∆ E thatcan be evaluated without explicitly converging the self-consistent procedure to the LS state. TABLE I: Magnetic exchange couplings (in meV) calculatedfrom energy derivatives and from energy differences for theH–He–H system. LSDA BLYP B3LYP J HS = S A S B “ d Edθ ” θ =0 − − − J LS = − S A S B “ d Edθ ” θ =180 ◦ − − − J ∆ E = E ( θ =180 ◦ ) − E ( θ =0)4 S A S B − − − Our second proof of concept was carried out in the ox-ovanadium(IV) dimer [( µ -OCH )VO(ma)] . This com-plex shows a strong antiferromagnetic coupling of about − Here we employedAhlrich’s triple-zeta valence basis set for for the Vatoms and Ahlrich’s double-zeta valence basis for first-row atoms , as obtained from Ref. 48. This basiswas shown to provide reliable results in practical cal-culations of exchange couplings.
Atomic coordinateswere taken from experimental crystallographic data. InFigs. 2 and 3 we present a plot of E KS as a functionof θ (0 ≤ θ ≤ ◦ ) for LSDA and B3LYP, respec-tively. In both figures, E KS ( θ ) follows closely a cosinefunction connecting the HS and LS extrema, indicatingthat both functionals capture the Heisenberg behaviorof the oxovanadium complex. Related investigations inperiodic systems using the LSDA can be found in theliterature. For the plots shown in Fig. 2 and Fig. 3 we have chosenas magnetic centers A and B ( S A = S B = 1 /
2) both setsof V and apical O atoms since most of the spin densityin this complex is localized on this moiety, as shown inFig. 4. However, it is worth to mention that by choosingthe metal atoms only as magnetic centers A and B thechanges in the plots are unappreciable and the calculated -45-40-35-30-25-20-15-10-5 0 0 20 40 60 80 100 120 140 160 180 E n e r gy ( m e V ) θ (degrees)High-Spin Low-Spinconstraint noncollinear spin DFT-2 J S S cos( θ ) FIG. 2: LSDA energy change as a function of the angle be-tween the local magnetic moments obtained from a constraintnoncollinear spin density functional calculation in the oxo-vanadium complex (Fig. 4). The solid line shows the (ideal)cosine function connecting the AF and FM extrema. -12-10-8-6-4-2 0 0 20 40 60 80 100 120 140 160 180 E n e r gy ( m e V ) θ (degrees)High-Spin Low-Spinconstraint noncollinear spin DFT-2 J S S cos( θ ) FIG. 3: Same as Fig. 2 for B3LYP. magnetic exchange couplings J HS and J LS vary very lit-tle. For instance, for LSDA, magnetic exchange couplingschange (in meV) from J HS = − . J HS = − . J LS = − . J LS = − . TABLE II: Magnetic exchange couplings (in meV) calcu-lated from energy derivatives and from energy differences forthe vanadium bimetallic complex. The corresponding exper-imental value is − a LSDA B3LYP J HS = S A S B “ d Edθ ” θ =0 − − J LS = − S A S B “ d Edθ ” θ =180 ◦ − − J ∆ E = E ( θ =180 ◦ ) − E ( θ =0)4 S A S B − − a Taken from Ref. 45 only 2.75 ◦ (for θ = 90 ◦ ). AB (a) AB (b) FIG. 4: Spin density isosurface of the HS (a) and LS (b)states of the oxovanadium complex (Fig. 4). Red correspondsto ↑ and blue to ↓ . The isosurface represents a spin densityof 0.01 (Bohr − ). Magnetic centers A and B chosen for thecalculations are also indicated. A careful comparison of Fig. 2 and Fig. 3 evidences alarger deviation from the ideal cosine function for LSDAthan for B3LYP. The largest differences from the cosinefunction are approximately 0.6 meV and 0.04 meV forLSDA and B3LYP, respectively, and occurs for θ = 90 ◦ in both cases. This is not surprising since LSDA yieldselectron (total and spin) densities more delocalized thanB3LYP (the L¨owdin atomic magnetic moments at the Vatoms for the LS state are 1.00 µ B and 1.10 µ B for LSDAand B3LYP, respectively) and therefore, one can expectthat the B3LYP energy follows the Heisenberg behaviormore closely than its LSDA counterpart. Localization ofthe spin-density also reduces the calculated magnetic ex-change couplings, as shown by Martin , Ruiz , anddemonstrated by Rudra et at. by explicitly constrainingthe local magnetization of the LS state. As shown inTable II, The difference between J HS and J LS is 4.7 meVfor LSDA, while it is only 0.2 meV for B3LYP. Thus, incontrast to the perfectly localized case, a more delocal-ized magnetization yields to larger deviations from theideal Heisenberg behavior and hence greater differencesbetween J HS and J LS and at the same time to largerexchange couplings. IV. CONCLUSIONS
We have proposed a method for the calculation of mag-netic exchange couplings from noncollinear spin densityfunctional calculations that employs the second deriva-tive of the electronic energy of a single state with respectto a parameter, Eqs. (6) and (7). Within this approachthere is no need to search for different self-consistent so-lutions of spin-states as it is commonly done in methodsbased on energy differences, such as the SP or NP meth-ods, Eqs. (2) and (3). Our method utilizes perturbationtheory for the evaluation of magnetic exchange couplingsand therefore, in combination with standard analytic sec-ond derivatives techniques, it can potentially be used to compute exchange couplings very efficiently, opening thepossibility of “black-boxifying” the extraction of mag-netic exchange couplings from density functional theorycalculations.Our proof of concept calculations show very promis-ing results. For the cases studied we found that ourmethod reproduces exchange couplings obtained fromthe spin-projected approach based on energy differences.As expected from physical grounds, the agreement be-tween both methods improves when the DFT descriptionof the interaction between the magnetic centers is moreHeisenberg-like. In this sense, the curve E KS ( θ ) providesa quantitative measure of how well the electronic systemmimics the behavior of an ideal Heisenberg model. V. ACKNOWLEDGMENTS
This research was supported in part by an awardfrom Research Corporation. J.E.P. acknowledges sup-port from the President’s Research Investment Fund(PRIF) and a start-up grant from Central Michigan Uni-versity. L. Noodleman, J. Chem. Phys. , 5737 (1981). E. Ruiz, P. Alemany, S. Alvarez, and J. Cano, J. Am.Chem. Soc. , 1297 (1997). J. Kortus, M. R. Pederson, T. Baruah, N. Bernstein, andC. S. Hellberg, Polyhedron , 1871 (2003). G. Christou, D. Gatteschi, D. N. Hendrickson, and R. Ses-soli, MRS Bull. , 66 (2000). P. Hohenberg and W. Kohn, Phys. Rev. B , 864 (1964). U. von Barth and L. Hedin, J. Phys. C: Solid State Phys. , 1629 (1972). J. Kortus, C. S. Hellberg, and M. R. Pederson, Phys. Rev.Lett. , 3400 (2001). J. Ribas-Arino, T. Baruah, and M. R. Pederson, J. Chem.Phys. , 044303 (2005). L. Noodleman and D. A. Case, Adv. Inorg. Chem. , 423(1992). L. Noodleman and E. R. Davidson, , 131 (1986). E. Ruiz, J. Cano, S. Alvarez, and P. Alemany, , 1391(1999). L. Noodleman, J. G. Norman, J. H. Osborne, A. Aizman,and D. A. Case, Journal of the American Chemical Society , 3418 (1985). E. Ruiz, A. Rodriguez-Fortea, J. Cano, S. Alvarez, andP. Alemany, , 982 (2003). M. Nishino, S. Yamanaka, Y. Yoshioka, and K. Yamaguchi,J. Phys. Chem. A , 705 (1997). I. Rudra, Q. Wu, and T. Van Voorhis, J. Chem. Phys. ,024103 (2006). I. Rudra, Q. Wu, and T. Van Voorhis, Inorg. Chem. ,10539 (2007). D. Dai and M.-H. Whangbo, J. Chem. Phys. , 2887(2001). A. E. Clark and E. R. Davidson, J. Chem. Phys. , 7382 (2001). E. R. Davidson and A. E. Clark, J. Phys. Chem. A ,7456 (2002). W. Kohn and L. J. Sham, Phys. Rev. A , 1133 (1965). A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov, andV. A. Gubanov, J. Magn. Magn. Mater. , 65 (1987). J. K¨ubler, K.-H. H¨ock, J. Sticht, and A. R. Williams, J.Phys. F: Met. Phys. , 469 (1988). L. Nordsr¨om and D. J. Singh, Phys. Rev. Lett. , 4420(1996). T. Oda, A. Pasquarello, and R. Car, Phys. Rev. Lett. ,3622 (1998). P. Kurtz, F. F¨orster, L. Nordsr¨om, G. Bihlmayer, andS. Bl¨ugel, Phys. Rev. B , 024415 (2004). S. Yamanaka, D. Yamaki, Y. Shigeta, H. Nagao, Y. Yosh-ioka, N. Suzuki, and K. Yamaguchi, Int. J. QuantumChem. , 664 (2000). P. H. Dederichs, S. Bl¨ugel, R. Zeller, and H. Akai, Phys.Rev. Lett. , 2512 (1984). J. Sticht, K.-H. H¨ock, and J. K¨ubler, J. Phys.: Condens.Matter , 8155 (1989). L. M. Sandratskii, Adv. Phys. , 47 (1998). J. E. Peralta, G. E. Scuseria, and M. J. Frisch, Phys. Rev.B , 125119 (2007). N. R. M. Mayer, S. Kruger, J. Chem. Phys. , 4411(2001). J. E. Peralta and G. E. Scuseria, J. Chem. Phys. , 5875(2004). F. Wang and T. Ziegler, J. Chem. Phys. , 12191 (2004). I. Malkin, O. L. Malkina, V. G. Malkin, and M. Kaupp, J.Chem. Phys. , 244103 (2005). M. K. Armbruster, F. Weigend, C. van W¨ullen, andW. Klopper, Phys. Chem. Chem. Phys. , 1748 (2008). P.-O. L¨owdin, J. Chem. Phys. , 365 (1950). E. Ruiz, J. Cirera, and S. Alvarez, Coord. Chem. Rev. ,2649 (2005). Gaussian Development Version, Revision F.02, M. J.Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A.Robb, J. R. Cheeseman, J. A. Montgomery, Jr., T. Vreven,G. Scalmani, B. Mennucci, V. Barone, G. A. Petersson, M.Caricato, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R.Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda,O. Kitao, H. Nakai, X. Li, H. P. Hratchian, J. E. Peralta,A. F. Izmaylov, K. N. Kudin, J. J. Heyd, E. Brothers, V.Staroverov, G. Zheng, R. Kobayashi, J. Normand, J. L.Sonnenberg, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega,J. C. Burant, J. M. Millam, M. Klene, J. E. Knox, J. B.Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts,R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C.Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski,S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K.Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J.V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski,B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Ko-maromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe,W. Chen, M. W. Wong, and J. A. Pople, Gaussian, Inc.,Wallingford CT, 2006. R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, J. Chem. Phys. , 650 (1980). S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys ,1200 (1980). A. Becke, Phys. Rev. A , 3098 (1988). C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B , 785(1988). A. D. Becke, J. Chem. Phys. , 5648 (1993). P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J.Frisch, J. Chem. Phys. , 11623 (1994), see also R. H.Hertwig and W. Koch, Chem. Phys. Lett.
345 (1997). Y. Sun, M. Melchior, D. A. Summers, R. C. Thompson,S. J. Rettig, and C. Orvig, Inorg. Chem. , 3119 (1998). A. Schafer, Horn, and R. Ahlrichs, J. Chem. Phys. ,2571 (1992). A. Schafer, C. Huber, and R. Ahlrichs, J. Chem. Phys. , 5829 (1994). D. J. Feller, J. Comp. Chem. , 1571 (1996). E. Ruiz, A. Rodriguez-Fortea, J. Tercero, T. Cauchy, andC. Massobrio, J. Chem. Phys. , 074102 (2005). P. Kurz, G. Bihlmayer, K. Hirai, and S. Bl¨ugel, Phys. Rev.Lett. , 1106 (2001). P. Novak, I. Chaplygin, G. Seifert, S. Gemming,and R. Laskowski, Comput. Mater. Sci. (2008),doi:10.1016/j.commatsci.2008.01.028. R. L. Martin and F. Illas, Phys. Rev. Lett.79