A general procedure to evaluate many-body spin operator amplitudes from periodic calculations: application to cuprates

This content has been downloaded from IOPscience. Please scroll down to see the full text. Download details: IP Address: 150.214.182.99 This content was downloaded on 27/02/2015 at 11:58 Please note that terms and conditions apply. A general procedure to evaluate many-body spin operator amplitudes from periodic calculations: application to cuprates View the table of contents for this issue, or go to the journal homepage for more 2007 New J. Phys. 9 369 (http://iopscience.iop.org/1367-2630/9/10/369) Home Search Collections Journals About Contact us My IOPscience The open–access journal for physics New Journal of Physics A general procedure to evaluate many-body spin operator amplitudes from periodic calculations: application to cuprates Ibério de P R Moreira1,4, Carmen J Calzado2, Jean-Paul Malrieu3 and Francesc Illas1 1Departament de Química Física and Institut de Química Teòrica i Computacional (IQTCUB), Universitat de Barcelona and Parc Científic de Barcelona, C/Martí i Franquès 1, E-08028 Barcelona, Spain 2Departamento de Química Física, Universidad de Sevilla, C/Prof. García González s/n, E-41012 Sevilla, Spain 3IRSAMC, Laboratoire de Physique Quantique, Université Paul Sabatier, 118 Route de Narbonne, F-31062 Toulouse-Cedex, France E-mail: [email protected], [email protected],[email protected] and [email protected] New Journal of Physics 9(2007) 369 Received 18 May 2007 Published 12 October 2007 Online at http://www.njp.org/ doi:10.1088/1367-2630/9/10/369 Abstract. A general procedure is presented which permits the form of an extended spin Hamiltonian to be established for a given magnetic solid and the magnitude of its terms to be evaluated from spin polarized, Hartree–Fock or density functional calculations carried out for periodic models. The computational strategy makes use of a general mapping between the energy of pertinent broken-symmetry solutions and the diagonal terms of the spin Hamiltonian in a local representation. From this mapping it is possible to determine not only the amplitude of the well-known two-body magnetic coupling constants between near-neighbor sites, but also the amplitudes of four-body cyclic exchange terms.Ascrutiny ofthe on-sitespin densities providesadditional information and control of the many broken-symmetry solutions which can be found. The procedure is applied to the La2CuO4, Sr2CuO2F2, Sr2CuO2Cl2and Ca2CuO2Cl2square lattices and the SrCu2O3ladder compound. It is shown that a proper description of the magnetic structure of these compounds requires that two- and four-body terms are explicitly included in the spin Hamiltonian. The implications for the interpretation of recent experiments are discussed. 4Author to whom any correspondence should be addressed. New Journal of Physics 9(2007) 369 PII: S1367-2630(07)50628-2 1367-2630/07/010369+25$30.00 © IOP Publishing Ltd and Deutsche Physikalische Gesellschaft 2 Contents 1. Introduction 2 2. Theoretical background and methodology 6 3. Computational details 13 4. Results 14 4.1. Spin Hamiltonian parameters for the 2D square lattices ............. 14 4.2. Spin Hamiltonian parameters for the SrCu2O3spin ladder compound ...... 17 4.3. Critical analysis of the results ........................... 19 5. Summary, conclusions and possible extensions 21 Acknowledgments 23 References 23 1. Introduction The discovery of the anomalous properties of high-Tcsuperconducting cuprates (HTCSs) in the late eighties has triggered a considerable interest in the crystal and electronic structure of these compounds from both experimental and theoretical points of view [1]–[5]. An enormous research effort on ceramic materials, mostly based on copper oxides, has been carried out continuously to try to improve the properties of known structures and synthetic pathways and has resulted in the synthesis of a wide variety of cuprates. The impressive richness of lowdimensional magnetic behavior of the different copper compounds can be, to a large extent, traced back to the stacking of the distorted CuX6octahedra (or CuX5pyramidal or CuX4planar units) in the lattice [4,6]. Most of the cuprate based materials are formed by almost independent CuO4units with distant apical ligands to complete the strongly distorted CuO4X2or CuO4X units and, hence, the structure is dominated by the link between CuO4units. Depending on the nature of the counter ions and on the number of links between the different CuO4units, different structures can be formed ranging from the typical lamellar two-dimensional (2D) structure of the HTCSs to many lower dimensional structures by different combinations of edgesharing and corner-sharing CuO4plaquettes (or CuO4units) that can give rise to spin ladders (e.g. Srn−1CunO2n−1series with n⩾2 [7]–[9]), zigzag spin chains (e.g. SrCuO2[10]) and quasi-1D systems (e.g. A2CuO2(A =Ca,Sr) [11,12] or Li2CuO2[13]) formed by edge-sharing CuO4units. The electronic ground state of this kind of materials is usually described by the open shell nature of the Cu2+ ions arranged in the CuO4units in which the Cu(3d9)atomic configuration gives rise to a dx2–y2type hole with the lobes pointing towards the O ions. The resulting Cu–O–Cu pathways range from ∼90◦to 180◦and they are responsible for the rich variety of low-dimensional magnetic structures dominated by moderate ferromagnetic (FM) to strong antiferromagnetic (AFM) interactions. From the theoretical point of view, these systems are strongly correlated in nature, making standard band theory techniques based on density functional theory (DFT) unable to accurately describe either their valence or low energy spectrum [14]. However, it has been shown that hybrid exchange-correlation functionals can provide reliable descriptions of strongly correlated transition metal magnetic systems ([15] and references therein). New Journal of Physics 9(2007) 369 (http://www.njp.org/) 3 Most of the HTCS materials have a lamellar structure in which strong AFM interactions take place along 180◦Cu–O–Cu bonds in edge sharing Cu4O4plaquettes leading to a 2D network of effective spin S=1/2 particles. These strong magnetic interactions observed in the HTCS are thought to be fundamental ingredients of the high-Tcsuperconductivity microscopic mechanism [5,16]. Since these compounds may be regarded as effective S=1/2 spin lattices, their low energy spectrum and collective properties are assumed to be governed by a simplified Heisenberg Hamiltonian as in equation (1) ˆ H=X hi,ji Jij ˆ Si·ˆ Sj−1 4,(1) which only accounts for the magnetic coupling Jij between nearest-neighbor (NN) centers i and j. This is in agreement with the widely accepted general picture for HTC superconductivity involving a ‘Heisenberg sea’ where holes or electrons are introduced by doping the perfect structures. However, it has been claimed that to fully understand the magnetic excitations, infrared and neutron scattering spectra of 2D [17]–[23] and spin ladder cuprates [24,25] it is necessary to extend the spin Hamiltonian as in equation (2) ˆ H=X i,j Jij ˆ Si·ˆ Sj−1 4 +X i,j,k,l Jring ˆ Si·ˆ Sjˆ Sk·ˆ Sl+ˆ Si·ˆ Slˆ Sj·ˆ Sk−ˆ Si·ˆ Skˆ Sj·ˆ Sl−1 16 +· · · .(2) In this expression the constants 1/4 and 1/16 have been introduced to define the zero of energy as that of the FM solution. In this spin model, the signs and amplitudes of the local intersite magnetic interactions govern the collective properties of a spin lattice. They appear as the basic ingredients of the effective spin Hamiltonian which in full generality involves not only the two-body exchange Jij but also other interactions, such as those represented by the four-body cyclic term Jring, or even higher-order terms. The largest two-body couplings are expected to occur between NN sites although next-NN (NNN) interactions may be non-negligible or even of the same order of magnitude (cf CuGeO3 system [26,27]). Regarding the four-body operator terms, they may be important in Cu4O4 plaquettes since their origin lies in the cyclic circulation of electrons around the ring. Their importance and that of analogous cyclic six-body effects have been pointed out in other types of half-filled band systems such as the πsystem of conjugated organic molecules [28]. Similar four-body terms are crucial to describe the ground state properties of 3He [29]. The direct determination of the amplitude of the many-body terms of an extended spin Hamiltonian such as the one in equation (2) from experiment is, in general, impossible. Spin ladders with a variety of intersite distances represent an especially difficult case. In general, the experimental determination of the coupling constants is based on a series of hypotheses about the negligibility of interactions between ‘remote’ sites. From these hypotheses, a given spin model is assumed and validated only from a numerical fit of the thermodynamic or spectral properties. It is customary to consider NN interactions only although this may be an excessive simplification and eventually can lead to contradictory estimates of the dominant couplings. This is precisely the case of the cuprate spin ladders for which conflicting values of the Jrung/Jleg ratios ranging from 0.5 to 1.0 were proposed. Interactions involving sites at a longer distance or New Journal of Physics 9(2007) 369 (http://www.njp.org/) 4 involving various sites had to be invoked to rationalize the different experimental results arising from different techniques. One may think for instance in inter-ladder interactions, diagonal terms in the plaquettes or four-body cyclic effects. Indeed, several theoretical studies seem to consistently indicate that four-body terms (Jring) are crucial [21]–[24]. Here, it is important to point out that very recently Toader et al [30] provided strong experimental evidence of the importance of the Jring term in La2CuO4with Jring ≈0.5J. This value is comparable to the pairing energies and strongly suggests that the resulting circulating currents could have an important role in the mechanism of superconductivity. This is in contrast with previous estimates of Jring for 2D and spin ladder cuprates, obtained either from indirect measurements or from numerical simulations with an extended Heisenberg model, which propose substantially smaller amplitudes with Jring ∼0.3J[19,21,22,24,31]. The origin of these discrepancies relies on the choice made by different authors for the magnitude of the other coupling constants in the spin Hamiltonian of equation (2). Hence, the Jring term as evaluated by Toader et al relies on NN and NNN coupling constants of J=111.8meV and Jd= −11.4meV, respectively, extracted from one of the various fittings of the magnon spectrum [22]. However, it is important to point out that this value for Jis smaller than another experimental estimate of 135±6meV obtained with a NN Heisenberg Hamiltonian [32]. Moreover, one should note that the present estimate of a FM Jddiffers from previous theoretical predictions [33] and from fitting to experimental measurements on materials with similar exchange paths [34]. However, one should also recognize that the theoretical study by Annett et al [33], takes only into account twoelectron processes to evaluate this Jdterm while, as clearly explained by Toader et al [30] and references therein, the FM character of Jdarises precisely from three-electron exchange processes around a plaquette, which are of the same order of magnitude as four-electron exchange processes around the same plaquette and which are not taken into account in [33,34]. The discussion above illustrates the difficulties faced by experimentalists when attempting to extract the magnitude of the important terms and the need for independent accurate and unbiased theoretical predictions. For the SrCu2O3ladder compound a similar situation is encountered; the recent Raman response experiments by Schmidt et al suggest Jrung =140meV, Jring/Jrung =0.2 and Jleg/Jrung =1.5 [35], in contrast with previous work indicating a more isotropic behavior between rungs and legs [36]. From the preceding discussion it is clear that an accurate prediction of the various magnetic interactions entering in the spin Hamiltonian as in equation (2) is not only highly desirable to understand the ground state properties of this kind of systems but urgently needed. One could, for instance, start from some approximate electronic Hamiltonian such as a single-band model involving only the magnetic centers or a two-band model involving also the electrons of the bridging ligands. However, this is likely to introduce even more problems since it is difficult to assess the accuracy of such a multiparametric approach. An alternative and straightforward way to an unbiased estimate is the direct evaluation of the amplitude of the relevant magnetic interactions on a realistic model system treating all the electrons in a large enough basis set. If practicable, this pragmatic approach will present the advantage of providing an independent, unbiased and consistent prediction of the amplitude of relevant exchange parameters that can solve the difficulties encountered by the usual fitting techniques used by experimentalists. Unfortunately, the direct calculation of the important terms is not a simple task since it also requires the use of a model for the solid. In a first approach, one may neglect translational symmetry and define a properly embedded finite cluster, a fragment of the periodic lattice with two, three or four magnetic sites with their coordination ligands and perform the best ab initio New Journal of Physics 9(2007) 369 (http://www.njp.org/) 5 (using the exact Hamiltonian) explicitly correlated calculations. This strategy, in particular when using the so-called difference dedicated configuration interaction (DDCI), has provided very consistent and reliable evaluations of the two-body interactions in a large series of perovskites, oxides and 2D cuprates [37]–[41]. However, this approach rapidly faces computational limits due to the need for very large embedded cluster models. This is especially the case when attempting to extract four-body terms in ladder compounds [42]–[45]. Alternatively, one may exploit the periodic symmetry of the crystal but in this case explicitly correlated calculations are not feasible and one must rely on spin-polarized mean-field type approaches. In such a case, only the FM solution may be considered as a valid approximation of an eigenstate of the problem with properly defined spin quantum numbers [46,48]. The solutions with lower values of the square of the total spin operator S2(lower magnetization) cannot be properly treated. Fixing a priori different localized periodic spin distributions one may obtain a set of distinct selfconsistent solutions of different energies. However, these solutions are not eigenfunctions of the S2operator; they are (spin) broken-symmetry solutions. We will show that assigning the energy of the broken-symmetry solutions to the expectation value of the corresponding spin distribution of a Heisenberg Hamiltonian enables one to obtain estimates of the magnetic coupling constants. Such approaches have been rather extensively used in the field of molecular magnetism and may also be employed on embedded clusters to extract the two-body terms [48,49]. The mean-field calculations may use the exact Hamiltonian, and the corresponding information comes either from unrestricted broken-symmetry Hartree–Fock (UHF) solutions or from the similar ones obtained by means of DFT approaches. In the first case, the resulting evaluations of the AFM couplings are severely underestimated [15], [50]–[52]. In the second one, the introduction of electron correlation effects improves the result but the numerical values are highly sensitive to the chosen exchange-correlation potentials. Notice that the local density approximation (LDA) and the different generalized gradient approximations (GGA) in DFT frequently lead to metallic solutions, closed shell in nature [15], which miss the main physical features of the low energy states and cannot provide any information about the magnetic couplings. Actually, the best results are obtained with hybrid functionals initially proposed by Becke [53,54], for which a percentage of Fock exchange is mixed with a given DFT exchange-correlation functional. However, the results are again strongly dependent on the amount of Fock exchange included in the exchange potential [55,56]. Previous calculations have shown that the best percentage, at least for NiO [15] and cuprates [57,58], is around 35% Fock. Here, it is worth pointing out that hybrid DFT calculations need to include the non-local Fock contribution which can be accurately computed when the total density is expressed in terms of Gaussian type orbitals (GTO) basis sets. Recently, however, it has been shown that a reasonable accuracy can be reached employing plane waves as well [59]. In the present work, we report a generalization of the periodic symmetry-broken (UHF or DFT) approach which permits one to obtain not only the various NN and NNN most important two-body terms, as initially done for MF2(M =Cu, Ni, Co, Fe, Mn) compounds [52,60], but to predict higher-order terms such as the four-body cyclic terms appearing in equation (2). To this end broken-symmetry hybrid DF computations are carried out for the La2CuO4, Sr2CuO2F2, Sr2CuO2Cl2and Ca2CuO2Cl22D square lattices and the SrCu2O3two-leg spin ladder. A general theoretical framework is presented that permits one to extract accurate values for the many-body terms of extended spin Hamiltonians from these periodic first-principle calculations. In addition, it will be shown that these periodic calculations confirm the importance of the FM interladder exchange and of the four-body intra-plaquette operator. New Journal of Physics 9(2007) 369 (http://www.njp.org/) 6 This paper is organized as follows: in section 2, we develop the relevant methodological aspects to extract the amplitude of the parameters of a general spin Hamiltonian from periodic all-electron calculations of the electronic structure of the systems. First, the logics leading to a mapping between spin eigenfunctions and symmetry-broken solutions is developed in a general scheme starting from the Hubbard model for a plaquette with few remarks on the calculated solutions. Section 3describes the computational details of the calculations. The results obtained for the planar cuprates and the ladder system are reported in section 4and, finally, in section 5, we present the conclusions of this investigation and further extensions of the theory. 2. Theoretical background and methodology Let us assume that one may calculate a mean field single determinant description of the FM state in a given unit cell. This wavefunction is usually expressed in terms of symmetry-adapted (delocalized) Bloch functions. An appropriate localizing transformation of the singly occupied functions will define atom-centered Wannier magnetic orbitals [61,62]. The FM solution for the system may be written as ΦSz,max =|core·abc · · · n|,(3) where core stands for the closed-shell part of the wavefunction and a, b, c, ...,nrepresent the magnetic orbitals localized on the different magnetic centers. In this picture there is one electron per site and one can refer to the corresponding determinant as a valence bond (VB) neutral form (see also the definition in [63]). In the following, the expectation energy hΦSz,max|ˆ H|ΦSz,maxi = E0will be taken as energy origin (E0=0). The nelectrons in norbitals define a half-filled band problem. Each of the possible distributions of the magnetic electrons in the magnetic orbitals defines a determinant and the set of these determinants spans an orthogonal VB basis set. The matrix elements of the exact Hamiltonian in this basis define a valence configuration interaction matrix. In order to understand the physical effects arising from the interaction between the αand β(or spin up and spin down) electrons, we shall consider first a finite set of neighbor sites, the periodicity being introduced in a subsequent step. One must first recall that the lowest energy determinants are the neutral ones (in the VB sense), for which any magnetic orbital is singly occupied. The ionic determinants set, presenting positively and negatively charged sites, are of higher energy. Let us consider a given VB neutral determinant as ΦI=core·a· · · h¯ i jk¯ lm · · · n,(4) built with the local magnetic orbitals of the FM solution and where the bars indicate, as usual, beta spin electrons (hence, iis synonymous of iαand ¯ iof iβ). The expectation energy (within the exact non relativistic Hamiltonian) of such a determinant is given by equation (5) hΦI|ˆ H|ΦIi − E0=X i∈α(I)X j∈β(I) Kij,(5) where α(I)and β(I)are the set of spin orbitals of αand βspin in ΦIand iαand jβthe corresponding spin orbitals. This expectation energy is higher than that ofΦSz,max since it misses the direct (positive) exchange integrals between the αand βspin orbitals. It is reasonable to assume that the nonzero Kij integrals concern NN and NNN site pairs only. These integrals New Journal of Physics 9(2007) 369 (http://www.njp.org/) 7 have been estimated in cuprates [42]–[45], [57,64,65] and play an important role in defining the final value of the magnetic coupling constant. However, for the perturbative development given below we consider that the difference among the different Kij integrals is small enough to be neglected in front of the excitation energies to the VB ionic determinant and a common zeroorder energy E0will be kept for the VB neutral determinants in the perturbative developments. The strategy will now consist of establishing a relation between the diagonal elements of the Heisenberg Hamiltonian and the energies of symmetry-broken UHF solutions. This is accomplished by means of suitable perturbation expansions. The Heisenberg Hamiltonian is obtained from the exact one by considering the space of all VB determinants ΦIas a model space [66] in the frame of the effective Hamiltonian theory of Bloch or des Cloizeaux [67,68]. Let us call ˆ P the projector on this subspace of dimension N ˆ P=X I=1,N |ΦIi hΦI|,(6) the effective Hamiltonian ˆ Heff is defined by the Neigen equations ˆ Heff|ˆ P9mi = Em|ˆ P9mi,m=1,N(7) where the vectors |9miare the Neigenvectors of ˆ H ˆ H|9mi=Em|9mi,m=1,N,(8) having the largest projections on the model space, and Emthe corresponding eigenvalues. Since all VB neutral determinants have the same spatial part, the effective Hamiltonian may be written in terms of spin operators and, in principle, it will involve many-spin operators. Due to the large energy gap between the neutral and VB ionic determinants the effective Hamiltonian may be approached through the quasi-degenerate perturbation theory (QDPT). Let us consider the diagonal matrix element associated with a determinant ΦI. To the second-order DΦI ˆ Heff ΦIE=DΦI ˆ H ΦIE+X α∈ionic DΦI ˆ HΦαEDΦα ˆ HΦIE E0−E0 α =E0+X iα, jβ Kij −2t2 ij U!,(9) where we have considered that the direct exchange integrals Kij are small in front of the energy differences between the neutral and the VB ionic configurations in the evaluation of the secondorder corrections. The determinant ΦIinteracts with all the determinants obtained by charge transfers between αand βsites leading to ionic determinants of the form Φi+j−I=a+ jaiΦI;DΦi+j−I ˆ H ΦIE=tij,(10) and Φi−j+I=a+ ¯ ia¯ jΦI;DΦi−j+I ˆ H ΦIE=tij.(11) The amplitude of the hopping integral tij rapidly decreases with the distance between sites iand j. Comparing hΦI|ˆ Heff|ΦIiwith equation (2) one obtains −Jij =2 Kij −2t2 ij U!,(12) New Journal of Physics 9(2007) 369 (http://www.njp.org/) 8 which is a rather well-known relation. Let us suppose now that one tries to optimize the energy of a symmetry-broken determinant Φ0 Iobtained from ΦIby relaxing the spin orbitals. To the first-order this optimization is equivalent to a CI between ΦIand the singly excited determinants leading to 9I=ΦI−X iα, jβ c(1) ij Φi+j−I+Φi−j+I=ΦI−X iα, jβ tij UΦi+j−I+Φi−j+I,(13) the spinorbital rotations i0=i−X jβ tij Uj,(14) j0=j−X iα tij Ui,(15) giving delocalization tails to the αspin orbital on site ion to the orbitals of the neighbor sites occupied by βspin electrons (or equivalently, c(1) ij =tij/U). From 9Ione can define a UHF single determinant Φ0 I=core·a0· · · h0¯ i0j0k0¯ l0m0· · · n0,(16) which only differs from 9Iby second-order effects Φ0 I≈9I+O(2)(17) and its energy only differs from the second-order corrected energy associated with 9Iby thirdorder terms DΦ0 I ˆ H Φ0 IE≈DΦI ˆ H ΦIE−X iα, jβ 2t2 ij U+O(3). (18) Hence, EUHF I0=hΦI|ˆ Heff |ΦIi+O(3)≈E0−1 2P iα, jβ Jij.(19) Equation (19) relates the calculated energy of the I’ broken symmetry solution for the exact Hamiltonian, relative to the FM state, to the corresponding energy expression of the Heisenberg Hamiltonian in equation (1) and it provides a first useful relationship to extract the relevant twobody magnetic coupling constants from the different broken-symmetry solutions. Since this relation is also valid when ΦIrepresents a periodic spin distribution, equation (19) provides a practical way to extract these effective parameters from periodic calculations. This is precisely the procedure followed in previous works [52,60]. However, it is possible to use the perturbation developments above to generalize the procedure in such a way that extraction of higher-order terms in the spin Hamiltonian as in equation (2) becomes straightforward. Before continuing the development a caveat is necessary since the Wannier functions as in equation (3) or (4) refer to the whole crystal. In practice, one has to limit the number of spin distributions to a finite number which is given by all possible spin permutations in a large enough supercell. Convergence of the results with respect to the model chosen can be indeed verified by extending the supercell. However, since the magnetic coupling constants are local parameters [37], the results are not expected to vary significantly upon enlarging the supercell used to extract the relevant parameters. Therefore, starting from nnon equivalent spin distributions ΦI,ΦJ, . . . , ΦM, one may, in principle, determine mquantities Jij and, of course, discriminate the leading terms from New Journal of Physics 9(2007) 369 (http://www.njp.org/) 15 Table 1. Energy expressions per Cu atom (relative to FM phase) for the magnetic solutions in figure 1used to extract J=Jij(NN),Jd=Jij(NNN)and Jring in the 2D La2CuO4, Sr2CuO2F2, Sr2CuO2Cl2and Ca2CuO2Cl2layered cuprates. The three entries for each magnetic phase show the corresponding calculated values (in meV) for the Fock-35, (B3LYP), and [UHF] potentials. Phase Energy expression La2CuO4Sr2CuO2F2Sr2CuO2Cl2Ca2CuO2Cl2 140.1 153.7 130.8 146.7 AFM J(182.8) (214.6) (182.0) (196.6) [31.0] [33.1] [26.2] [32.4] 78.9 84.7 70.1 79.1 AFM1 J/2+ Jd/2+ Jring/8 (102.5) (125.1) (103.1) (109.4) [15.4] [16.9] [13.2] [16.3] 78.8 82.2 68.2 78.3 AFM2 J/2+ Jd(100.3) (118.1) (99.4) (105.5) [15.1] [17.0] [13.3] [16.0] 0.0 0.0 0.0 0.0 FM 0 (0.0) (0.0) (0.0) (0.0) [0.0] [0.0] [0.0] [0.0] coupling amplitudes (J) it is worth noting that UHF largely underestimates them, as found in many previous works [57,76,80]. Density functional calculated values generally represent an improvement over UHF results although they exhibit the typical strong dependence on the exchange-correlation functional already pointed out by Martin and Illas [55,56]. Hence, B3LYP overestimates the Jvalues by a factor close to 1.5 whereas the values predicted by the Fock-35 potential are close to experiment as also expected from previous work on similar systems [15, 76,80]. Therefore, one can take the Fock-35 values as a reliable prediction. In all cases, the NNN diagonal interaction is predicted to be small and of AFM character, in contrast to what is suggested by some fits to experimental data [81]. However, the most important result of the present work concerns the cyclic exchange Jring amplitude. The general procedure described in detail in section 2has permitted the first direct estimate of the amplitude of the Jring terms from ab initio periodic calculations [82]. For La2CuO4, the outcome of the present periodic approach using the Fock-35 parameterization is in agreement with the other only available theoretical estimates for this compound arising from cluster calculations performed by either explicitly correlated wavefunctions [42,43] or DFT based calculations [45,76] and experimental estimations of J[22,32]. In addition, the present Fock-35 estimate of Jring ∼35.8meV is in excellent agreement with the indirect evaluations of Coldea et al (Jring ∼38±8meV) [22] and the simulations of Mizuno et al (Jring ∼40meV) [31]. Notice, however, that the Jring/J=0.25 ratio predicted here is smaller than that reported by Toader et al [30], which is of 0.5, but in excellent agreement with a more generally accepted ratio Jring/J=0.3 [21,22,31] for this compound. Finally, we would like to point out that the experimental determination of this term is still under discussion and even the existence and importance of four-body terms in La2CuO4have raised some controversy [83]. The present work provides unbiased first principles results that fully New Journal of Physics 9(2007) 369 (http://www.njp.org/) 16 Table 2. Numerical values (in meV) for the two- and four-body exchange amplitudes in La2CuO4, Sr2CuO2F2, Sr2CuO2Cl2and Ca2CuO2Cl22D cuprates obtained using the Fock-35, (B3LYP), and [UHF] potentials. Available experimental data are discussed in sections 1and 4.1. System J JdJring Jring/J Tc(max)(K) 140.1 8.8 35.8 0.25 La2CuO4(182.8) (8.9) (53.1) (0.29) 42a [30.9] [−2.8] [0.7] [0.02] 153.7 5.4 41.1 0.27 Sr2CuO2F2(214.6) (10.7) (99.5) (0.46) 46b [33.1] [0.4] [1.4] [0.04] 130.0 2.8 26.4 0.20 Sr2CuO2Cl2(182.0) (8.4) (62.8) (0.34) – [26.2] [0.1] [0.2] [0.01] 146.7 4.9 26.3 0.18 Ca2CuO2Cl2(196.6) (7.2) (60.0) (0.31) 28c [32.4] [−0.2] [1.4] [0.04] aFrom [87], b[71], c[88]. support the arguments of Toader et al in their reply to Raymond et al [84]. It is also worth pointing out that the Fock-35 calculations predicts an AFM Jd∼8.8meV again consistent with the values provided by embedded cluster calculations [42,44,77]. Hence, it will be of great interest to repeat the fit in Toader et al experiments [30] by using the present estimate of both Jand of Jdor, in a bottom up approach, use the present values for the three effective parameter to check consistency with respect to experiment. For Ca2CuO2Cl2, Sr2CuO2F2and Sr2CuO2Cl2, there are not previous theoretical or experimental values for the amplitude of the ring exchange term. Present calculations predict values for Ca2CuO2Cl2and Sr2CuO2Cl2that are slightly smaller than the corresponding ones in La2CuO4. This is consistent with larger NN and NNN distances in Ca2CuO2Cl2and Sr2CuO2Cl2 compared to La2CuO4(3.87 and 3.97 versus 3.81Å, for NN and ∼5.47 and ∼5.62 versus 5.38Å for NNN distances). For Sr2CuO2Cl2the agreement between the estimate of the NN interaction (J∼130meV) and the experimental value (J=125±6meV [85,86]) is excellent, the Jdcoupling is also predicted to be AFM and the estimated Jring/Jratio ∼0.20 is somewhat smaller than that for La2CuO4but still in a realistic range. For Sr2CuO2F2the results are closer to La2CuO4ones despite the differences in crystal structure (3.86 versus 3.81Å, for NN and ∼5.45 and ∼5.62 versus 5.38Å for NNN distances). Before closing this section it is interesting to note that the magnitude of Jalone does not discriminate Sr2CuO2F2and Ca2CuO2Cl2from La2CuO4(stoichiometric Sr2CuO2Cl2 corresponds to a extremely stable structure and attempts to synthesize doped phases by cation substitution of interstitial anion excess have been unsuccessful). Taking the reliable Fock-35 values, for La2CuO4Tc=42K [87] and J=140.1meV and, similarly, Sr2CuO2F2Tc= 46K [71] but J=153.7meV whereas for Ca2CuO2Cl2Tc=28K [88] but J=146.7meV. A larger Jvalue does not correspond to a larger Tc. However, Jring points towards a different New Journal of Physics 9(2007) 369 (http://www.njp.org/) 17 behavior since Jring =35.8meV for La2CuO4and Jring =41.1meV for Sr2CuO2F2whereas Jring =26.3meV for Ca2CuO2Cl2. This is in line with the results of Toader et al suggesting that the Jring plays an important role in defining the properties of the superconducting phase. We would like to point out that this is also in line with previous results indicating that a linear relationship exists between J/tand Tc[41,76]. This is because, a correct value of J/talso indicates a correct value of t/Uand hence of Jring/J(∼(t/U)2, cf equation (22) and [45]). 4.2. Spin Hamiltonian parameters for the SrCu2O3spin ladder compound For the SrCu2O3spin ladder, ten different magnetic solutions have been obtained from the spin distributions shown in figure 2. Due to an excessive spin frustration between the ladders, the FM solution becomes exceedingly high in energy and indeed physically meaningless. The shortrange intra-ladder repulsion forces the unpaired electron to be excessively delocalized on the oxygen atoms. Therefore, the AFM gs solution is used as energy origin. For the remaining nine broken-symmetry solutions with different magnetic orders, we derived a set of nine equations which are shown in table 3. These equations have been derived assuming that only five magnetic interactions may have non-negligible amplitudes: these are NN interactions along the legs Jl, or the rungs Jr, the NNN Jdinteractions in the plaquettes, the interaction between legs of different ladders Ji, and the four-body operator Jring amplitude. At variance from the abovediscussed cuprates, the resulting equations are linearly dependent, the system of equations is overdetermined and a least-square procedure has been used to find the optimum set of values for the five magnetic interactions considered. The optimum values are reported in table 4. The consistency of the set of equations is almost complete since the standard deviation between the DF calculated energies and those obtained from the equations in table 3and the so-obtained set of magnetic interactions is of 1.5meV for the B3LYP results and even smaller for the Fock-35 set. This consistency indicates that all non-negligible interactions have been included in the spin model Hamiltonian. For all the methods the dominant interactions are the couplings across the rungs (Jr) and along the legs (Jl) which all methods predict to be nearly equal in magnitude—i.e. (Jr/Jl)≈1— as expected from the similarity of the Cu–O–Cu exchange paths in agreement with previous two magnetic sites cluster model calculations [36] although a moderate cluster size dependence on the magnetic interaction along the leg has been observed [77]. Nevertheless, enlarging the cluster model by introducing a third copper atom belonging to the very close neighbor ladder results again in Jr/Jl≈1 [71]. In any case, this dependence evidences the difficulty of embedding finite clusters for certain types of materials and supports the present periodic approach that does not depend on the cluster design used to represent the real material. The accuracy of the present calculations can be judged from the amplitude predicted for Jrthrough the most realistic Fock-35 potential (∼155meV) which is consistent with the recent Raman response experiments [35]. The estimated Jdvalue is AFM as in the 2D cuprates and a non-negligible interladder FM exchange (Ji≈ −0.22Jl) is found and, consequently, one should not consider this system as formed by non-interacting ladders. Regarding Jring, both its amplitude and the Jring/Jrratio are larger than for 2D cuprates (table 2), but still close to 0.3 in agreement with various sets of experimental informations [21]–[24]. The agreement is particularly good with the values reported by Schmidt et al [35] although these authors do not consider Ji. This neglect is probably the reason for the strong anisotropy in their fitting which results in a Jl/Jr=1.5 ratio, a value which is difficult to accept from the Cu–O–Cu distances. New Journal of Physics 9(2007) 369 (http://www.njp.org/) 18 Figure 2. Schematic representation of the magnetic solutions calculated to extract Jl,Jr,Jd,Jiand Jring parameters in the ladder cuprate SrCu2O3. Continuous (dashed) lines correspond to FM (AFM) alignments for Jland Jr in the ladders. The gray zone defines the cell used for the symmetry-broken solution. Such large anisotropy between Jland Jris not supported from the present periodic calculations. We also found that repeating the least square fitting of the energies in table 2but neglecting the terms involving Jring results in a larger anisotropy in the thus Jrand Jlcalculated values. These results strongly suggest revising previous fittings using the present estimates as starting points. New Journal of Physics 9(2007) 369 (http://www.njp.org/) 19 Table 3. Energy expressions per Cu atom (relative to AFM gs phase) for the magnetic solutions calculated to extract Jl,Jr,Jd,Jiand Jring in the SrCu2O3 ladder compound (for the definition of these interactions see section 2). These correspond to the magnetic solutions given schematically in figure 2. The last column shows the corresponding calculated energy values (in meV) for the Fock-35, (B3LYP), and [UHF] potentials. Phase Energy expression Values in SrCu2O3 AFM9 Jl/6+ Jr/6−3Jd/2−0.05Jring 44.98 (57.15) [10.79] AFM8 Jl/4−Jd/4+ Ji/8 31.92 (48.85) [6.86] AFM7 Jl/2−Jd/2+ Ji/4 67.38 (96.23) [14.23] AFM6 Jl/2−Jd/2−Ji/4 83.78 (114.37) [18.53] AFM5 Jl/4+ Jr/8−Jd/4−Ji/8−0.075Jring 57.68 (75.75) [13.35] AFM4 Jr/4−Jd/2 38.87 (50.22) [8.14] AFM3 Jl/4+ Jr/8−Jd/4−0.075Jring 53.43 (70.73) [12.24] AFM2 Jr/8−Jd/4 16.20 (20.41) [4.07] AFM1 Jl/8+ Jr/16−Jd/8−0.0375Jring 26.68 (35.39) [6.15] AFM gs 0 0.0 (0.0) [0.0] Table 4. Numerical values (in meV) for the two- and four-body exchange amplitudes in SrCu2O3ladder compound obtained using the Fock-35, (B3LYP), and [UHF] potentials. Available experimental data are discussed in sections 1 and 4. System JlJrJdJiJring Jring/Jr 153.1 155.6 2.7 −34.2 48.8 0.31 SrCu2O3(216.3) (204.1) (5.4) (−36.0)(97.0) (0.47) [32.7] [32.7] [0.1] [−8.8] [0.04] [ < 0.01] 4.3. Critical analysis of the results Tables 2and 4summarize the results for the effective spin Hamiltonian parameters of the 2D cuprates and the ladder compound, respectively, obtained by means of the UHF method, B3LYP and Fock-35 hybrid potentials. Since the Fock-35 values are in agreement with available experimental data it is clear that B3LYP values are overestimations. This result merits a further deeper analysis. As well known, the NN magnetic coupling Jscales as t2/U(cf equation (12)). Previous works have shown that tvalues are not very sensitive to the exchange potential [57,77]. This suggests that the exceedingly large Jvalue predicted by the B3LYP method effectively implies a too small Uvalue and hence an insufficient on-site two-electron repulsion, as discussed elsewhere [57,80]. This incorrect description of electron–electron correlations in the B3LYP functional has a much more dramatic consequence. In fact, it has been shown by Malrieu and Maynau [28] that Jring scales as t4/U3(cf equation (22)) and, hence, a too small U value will result in an even larger overestimation of Jring as clearly shown in tables 2and 4. This interpretation is supported by the results obtained from the UHF method. In this case the values New Journal of Physics 9(2007) 369 (http://www.njp.org/) 20 La2CuO4Ca2CuO2Cl2 Sr2CuO2F2 Sr2CuO2Cl2 Figure 3. Representation of atomic spin densities on center i(ρi(I)), defined according to Mulliken population analysis, versus the number of opposed spins on NN sites (nNNβ(I)) of the different Imagnetic ‘phases’ for the 2D cuprates using UHF, Fock-35 and B3LYP methods. A linear fit, as suggested by equation (21), is also shown. of Jare largely underestimated due to the lack of electronic correlation effects, a well-known effect [39,55,56,90,91]. Therefore, one may conclude that UHF largely overestimates the on site two-electron repulsion and, consequently, Jring obtained at this level of theory must be very small. This is indeed the case, the UHF estimate of Jring is vanishingly small. These results are manifestations of the sensitivity of Jring/Jwith respect to the |t/U|ratio; this is Jring/J scaling as |t/U|2. If J(calc)/J(exp)=λ, then Jring(calc)/Jring(exp)≈λ3. To conclude, the DFT based approaches effectively include electronic correlation effects that lead to a decrease of the effective Uparameter which was severely overestimated by the mean-field Hartree–Fock method. However, the B3LYP potential appears to lead to unphysical too low Uvalues with important consequences in the resulting physical description. Scrutiny of the atomic spin densities, defined according to Mulliken population analysis, brings consistent additional information. In order to check the validity of equation (21), we report in figure 3the dependence of the spin densities on the metal centers as a function of the number of opposed spins on NN sites. Notice that this information is extracted from the New Journal of Physics 9(2007) 369 (http://www.njp.org/) 21 various magnetic ‘phases’ for each compound. As suggested by equation (20), a good linear dependence of the local spin density on the number of spin alternations is observed for each of the Hamiltonians. This result evidences the different strengths of the electronic delocalization in these three Hamiltonians and the slope of the correlation can be used to estimate the t/U ratio. For La2CuO4, the t/Uestimated ratios are 0.044, 0.093 and 0.110 for UHF, Fock-35 and B3LYP, respectively. These results nicely compare with previous estimates of this ratio from finite cluster calculations: 0.080 and 0.120 for Fock-35 and B3LYP [57]. Also, these results confirm the general trend discussed above, namely an excessive localization in UHF, resulting in a too large unscreened effective U, and a slightly too strong delocalization in B3LYP with a concomitant underestimation of the effective U. In the ladder compound, the analysis of the correlation with respect to the number of spin alternations is not univocal since the bonds are not identical and, therefore, is not carried out. As mentioned in section 2, it is observed that some highly frustrated and high energy selfconsistent solutions provide local spin densities which deviate significantly from the behaviour expected from equation (21). This is due to the physically meaningless strong metal-to-ligand delocalization of the magnetic electrons. Hence, the spin density analysis provides a criterion to eliminate these spurious solutions. This problem has been essentially observed in the ladder where the Jvalues are larger as a result of the shorter Cu–Cu distances and, also, due to the FM interladder interaction. As a consequence, the total spin frustration becomes exceedingly costly. Accordingly, we strongly recommend using a sufficiently large set of low energy solutions in the fitting procedure. 5. Summary, conclusions and possible extensions The extraction of NN spin couplings in magnetic lattices from symmetry-broken mean-field calculations has been frequently performed on molecular systems and solids (see [48] and references therein). The present work uses the benefit of the large multiplicity of symmetrybroken periodic solutions and proposes a new, general and unbiased scheme to predict the amplitude of the parameters defining a general spin Hamiltonian from DFT periodic calculations. An important property of the present procedure is that one does not need to make any assumption on the relative amplitude of these terms. Instead, it relies on a mapping approach between the energy of pertinent magnetic solutions and the diagonal terms of the spin Hamiltonian in a local representation. In addition, the overdetermined set of equations provides a test for the completeness of the interactions considered in the spin model. It has also been shown that the spin densities of the symmetry-broken solutions may be rationalized and they furnish additional information on the t/Uratio through equation (20). Hence, the analysis of the spin density distributions of these various solutions provides useful and consistent information on the t/U(or J/t) ratio and, hence, on the electronic delocalization. As a corollary one can point out that the observed correlation between Tcat optimal doping and the J/tratio [41, 76] can be formulated as well as a correlation with t/Uratio that brings in a more physically intuitive description. Finally, this work also points out the extreme sensitivity of the four-body operators amplitudes to the choice of the density functional and especially to the ratio of Fock exchange. For the time being the use of a 35% ratio has led to consistent results for cuprates and other transition metal compounds. The question of the transferability of this parametric quantity to other magnetic ions remains open. New Journal of Physics 9(2007) 369 (http://www.njp.org/) 22 A very important outcome of the present work concerns the prediction of the four-body cyclic operators in the general spin Hamiltonian in equation (2). The application of the general procedure discussed at length in section 2to a series of compounds provides an independent confirmation of the importance of four-body terms in these materials. In particular, for La2CuO4 the present results provide further support to the conclusions by Toader et al [30] based on the analysis of neutron diffraction experiments. Moreover, we supply reliable values for spin Hamiltonian parameters of other key compounds such as Ca2CuO2Cl2, Sr2CuO2F2and Sr2CuO2Cl2. These permit us to conclude that the importance of four-body terms is likely to be similar for most of the HTCS related cuprates and even more important in ladder compounds, in agreement with the recent experiments of Schmidt et al [35]. Likewise, it is also suggested that the fit to neutron scattering data should be revised by considering alternative values for both Jand Jdmagnetic coupling terms. In addition, the present study provides further evidence that the four-body term and the FM inter-ladder exchange introducing spin frustration between legs in the SrCu2O3ladder should not be neglected. To conclude this study, it is worth mentioning some possible generalizations of the procedure presented in section 2. One may, for instance, wonder whether the six-body spin operators would have significant amplitudes in such lattices. The cyclic circulation of electrons in six-member rings—related to the so-called chemical aromaticity—is responsible for the appearance of such operators at the sixth-order of perturbation theory. It has been demonstrated [28] that while the four-body operator results in a coupling = 40 t 4 /U 3 H ˆ the six-body matrix element resulting in a rotation of six spins is given by = 504 t 6 /U . 5 H ˆ The large prefactor is due to the huge number of processes leading to a full permutation of spins in a six-membered ring. Now, notice that a set of two fused plaquettes defines such a sixmembered ring. It would be worth checking whether, despite the smallness of the t/Uratio, such many-body operators are negligible in cuprate lattices. From the perturbative arguments above, and in view of the calculated t/Uratio, one expects the ratio Jsix-ring/Jfour-ring ∼0.1. Consequently, this question has not been addressed in the present study. Yet another generalization would concern lattices with S>1/2 magnetic sites such as Ni(d8)ions, NiO being the paradigm of these materials. The same strategy is applicable to evaluate not only the NN interactions but four-body exchange amplitudes as well. This will be considered in further work which will also show that one may, in principle, use symmetry-broken solutions to evaluate the possible importance of biquadratic spin exchange terms. New Journal of Physics 9(2007) 369 (http://www.njp.org/) 23 Acknowledgments Financial support from the French CNRS (UMR5626), Spanish Ministerio de Educacion y Ciencia (Ramon y Cajal research contract (I de PRM) and projects CTQ2005-08459-CO2-01, UNBA05-33-001), the Generalitat de Catalunya (projects 2005SGR00697, 2005 PEIR 0051/69 and Distinció per a la Promoció de la Recerca Universitària granted to FI), European Community (project COST-D26/0008/02) and French–Catalan international cooperation project (PICS 1458) is fully acknowledged. References [1] Mattis D C 2006 The Theory of Magnetism Made Simple (Singapore: World Scientific) [2] Kahn O 1993 Molecular Magnetism (New York: Wiley) [3] Yosida K 1998 Theory of Magnetism, Springer Series in Solid-State Sciences vol 122 (Berlin: Springer) [4] Cava R J 2000 J. Am. Ceram. Soc. 83 5 [5] Dagotto E 1994 Rev. Mod. Phys. 66 763 [6] Rice T M 1998 Physica B 241 [7] Hiroi Z, Azuma M, Takano M and Bando Y 1991 J. Solid State Chem. 95 230 [8] Takahashi H, Mori N, Nakanishi T, Nagata T, Uehara M, Akimitsu J and Kinoshita K 1997 Physica B 112 237 [9] Azuma M, Hiroi Z, Takano M, Ishida K and Kitaoka Y 1994 Phys. Rev. Lett. 73 3463 [10] Teske C L and Müller-Buschbaum H 1970 Z. Anorg. Allg. Chem. 379 234 [11] Hjorth M and Hyldtoft J 1990 Acta Chem. Scand. 44 516 [12] Ami T, Crawford M K, Harlow R L, Wang Z R, Johnston D C, Huang Q and Erwin R W 1995 Phys. Rev. B 51 5994 [13] Sapiña F, Rodríguez-Carvajal J, Sanchís M J, Ibañez R, Beltrán A and Beltrán D 1990 Solid State Commun. 74 779 [14] Pickett W E 1989 Rev. Mod. Phys. 61 433 Pickett W E 1989 Rev. Mod. Phys. 61 749 (correction) [15] Moreira I de P R, Illas F and Martin R L 2002 Phys. Rev. B65 155102 [16] Nagaosa N 1997 Science 275 1078 [17] Takahashi M 1977 J. Phys. C: Solid State Phys. 10 1289 [18] Roger M and Delrieu J M 1989 Phys. Rev. B39 2299 [19] Schmidt H and Kuramoto Y 1990 Physica C167 263 [20] Honda Y, Kuramoto Y and Watanabe T 1993 Phys. Rev. B47 11329 [21] Lorenzana J, Eroles J and Sorella S 1999 Phys. Rev. Lett. 83 5122 [22] Coldea R, Hayden S M, Aeppli G, Perring T G, Frost C D, Mason T E, Cheong S-W and Fisk Z 2001 Phys. Rev. Lett. 86 5377 [23] Sugai S, Sato M, Kobayashi T, Akimitsu J, Ito T, Takagi H, Uchida S, Hosoya S, Kajitani T and Fukuda T 1990 Phys. Rev. B42 R1045 [24] Brehmer S, Mikeska H-J, Müller M, Nagaosa N and Uchida S 1999 Phys. Rev. B60 329 [25] Matsuda M, Katsumata K, Eccleston R S, Brehmer S and Mikeska H-J 2000 Phys. Rev. B62 8903 [26] Castilla G, Chakravarty S and Emery V J 1995 Phys. Rev. Lett. 75 1823 [27] Fabricius K, Klümper A, Löw U, Büchner B, Lorenz T, Dhalenne G and Revcolevschi A 1998 Phys. Rev. B 57 1102 [28] Malrieu J P and Maynau D 1982 J. Am. Chem. Soc. 104 3021 [29] Roger M, Delrieu J M and Hetherington J H 1983 Rev. Mod. Phys. 55 1 [30] Toader A M, Goff J P, Roger M, Shannon N, Stewart J R and Enderle M 2005 Phys. Rev. Lett. 94 197202 New Journal of Physics 9(2007) 369 (http://www.njp.org/) 24 [31] Mizuno Y, Tohyama T and Maekawa S 1998 Phys. Rev. B58 R14713 Mizuno Y, Tohyama T and Maekawa S 1999 J. Low Temp. Phys. 117 389 [32] Kreimer B et al 1992 Phys. Rev. B46 14034 [33] Annett J F, Martin R M, McMahan A K and Satpathy S 1989 Phys. Rev. B40 2620 [34] Kim Y J et al 1999 Phys. Rev. Lett. 83 852 [35] Schmidt K P, Gössling A, Kuhlmann U, Thomsen C, Löffert A, Gross C and Assmus W 2005 Phys. Rev. B 72 094419 [36] de Graaf C, Moreira I de P R, Illas F and Martin R L 1999 Phys. Rev. B60 3457 [37] Moreira I de P R, Illas F, Calzado C J, Sanz J F, Malrieu J P, Ben Amor N and Maynau D 1999 Phys. Rev. B 59 R6593 [38] Calzado C J, Sanz J F, Malrieu J P and Illas F 1999 Chem. Phys. Lett. 307 102 [39] Calzado C J, Sanz J F and Malrieu J P 2000 J. Chem. Phys. 112 5158 [40] Muñoz D, Illas F and Moreira I de P R 2000 Phys. Rev. Lett. 84 1579 [41] Moreira I de P R, Muñoz D, Illas F, de Graaf C and Garcia-Bach M A 2001 Chem. Phys. Lett. 345 183 [42] Calzado C J and Malrieu J P 2001 Phys. Rev. B63 214520 [43] Calzado C J and Malrieu J P 2001 Eur. Phys. J. B21 375 [44] Calzado C J, de Graaf C, Bordas E, Caballol R and Malrieu J P 2003 Phys. Rev. B67 132409 [45] Calzado C J and Malrieu J P 2004 Phys. Rev. B69 094435 [46] Illas F, Moreira I de P R, Bofill J M and Filatov M 2004 Phys. Rev. B70 132414 [47] Illas F, Moreira I de P R, Bofill J M and Filatov M 2006 Theor. Chem. Acc. 115 587 [48] Moreira I de P R and Illas F 2006 Phys. Chem. Chem. Phys. 81645 [49] Moreira I de P R and Illas F 1997 Phys. Rev. B55 4129 [50] Ricart J M, Dovesi R, Roetti C and Saunders V R 1995 Phys. Rev. B52 2381 Ricart J M, Dovesi R, Roetti C and Saunders V R 1997 Phys. Rev. B55 15942 (erratum) [51] Dovesi R, Ricart J M, Saunders V R and Orlando R 1995 J. Phys.: Condens. Matter 77997 [52] Reinhardt P, Moreira I de P R, de Graaf C, Illas F and Dovesi R 2000 Chem. Phys. Lett. 319 625 [53] Becke A D 1993 J. Chem. Phys. 98 1372 [54] Becke A D 1988 Phys. Rev. A38 3098 [55] Martin R L and Illas F 1997 Phys. Rev. Lett. 79 1539 [56] Illas F and Martin R L 1998 J. Chem. Phys. 108 2519 [57] Calzado C J, Cabrero J, Malrieu J P and Caballol R 2002 J. Chem. Phys. 116 3985 [58] Moreira I de P R and Dovesi R 2004 Int. J. Quantum Chem. 99 811 [59] Franchini C, Bayer V, Podloucky R, Paier J and Kresse G 2005 Phys. Rev. B72 045132 [60] Moreira I de P R, Dovesi R, Roetti C, Saunders V R and Orlando R 2000 Phys. Rev. B62 7816 [61] Wannier G H 1937 Phys. Rev. 52 191 [62] Harrison W A 1989 Electronic Structure and the Properties of Solids (New York: Dover) p 87 [63] Moreira I de P R, Suaud N, Guihéry N, Malrieu J P, Caballol R, Bofill J M and Illas F 2002 Phys. Rev. B66 134430 [64] Stechel E B and Jennison D R 1988 Phys. Rev. B38 4632 [65] Hybertsen M S, Stechel E B, Schluter M and Jennison D R 1990 Phys. Rev. B41 11068 [66] Durand P, Malrieu J P and Lawley K P 1987 Ab Initio Methods in Quantum Chemistry vol 1 (New York: Wiley) chapter 4 [67] Bloch F 1958 Nucl. Phys. 6329 [68] des Cloizeaux J 1960 Nucl. Phys. 20 321 [69] Miller L L, Wang X L, Wang S X, Stassis C, Johnston D C, Faber J Jr and Loong C K 1990 Phys. Rev. B41 1921 [70] Argyriou D N, Jorgensen J D, Hitterman R L, Hiroi Z, Kobayashi N and Takano M 1995 Phys. Rev. B51 8434 [71] Al-Mamouri M, Edwards P P, Greaves C and Slaski M 1994 Nature 369 382 New Journal of Physics 9(2007) 369 (http://www.njp.org/)