Magnetic and Electronic Properties of Complex Oxides from First‐Principles

The theoretical treatment of complex oxide structures requires a combination of efficient methods to calculate structural, electronic, and magnetic properties, due to special challenges such as strong correlations and disorder. In terms of a multicode approach, this study combines various complementary first‐principles methods based on density functional theory to exploit their specific strengths. Pseudopotential methods, known for giving reliable forces and total energies, are used for structural optimization. The optimized structure serves as input for the Green's function and linear muffin‐tin orbital methods. Those methods are powerful for the calculation of magnetic ground states and spectroscopic properties. Within the multicode approach, disorder is investigated by means of the coherent potential approximation within a Green's function method or by construction of special quasirandom structures in the framework of the pseudopotential methods. Magnetic ground states and phase transitions are studied using an effective Heisenberg model treated in terms of a Monte Carlo method, where the magnetic exchange parameters are calculated from first‐principles. The performance of the multicode approach is demonstrated with different examples, including defect formation, strained films, and surface properties.

therefore attracted significant attention. Of special interest are multiferroic materials, [9] because several ferroic phases are coupled nontrivially giving rise to novel tuning mechanisms and control of future electronic and spintronic devices. Another significant breakthrough in oxide materials research was achieved by the discovery of the first oxide quasicrystal (OQC), [10,11] where it was shown that bariumtitanate, intensively investigated and applied over decades, shows a quasicrystalline structure as a monolayer on a Pt (111) surface. This underlines the importance of the study of oxide thin films.
Oxides containing transition metals or rare earth elements furthermore exhibit strongly localized electrons. Thus, the weak overlap of localized orbitals with orbitals of neighboring atoms leads to a decrease of the electron hopping strength resulting in flat bands. In such flat bands the involved energy range in terms of the band width becomes comparable with the electronelectron interaction strength leading to strong correlation effects and results in rich phase diagrams containing various phases of magnetic, charge, or orbital ordering as well as unconventional superconductivity. [12,13] The general development of novel oxide-based devices and oxide materials research-only a few trends were mentioned earlier-revealed several challenges in basic research. [14] In addition to understanding the underlying physics of the emerging effects, a major bottleneck remains in identifying respective functional materials for future technical applications. This concerns the identification of relevant oxides, the growth and fabrication of oxide thin films, and the development of control and tuning mechanisms based on strain and defect concentration. A combination of first-principles methods based on density functional theory (DFT) and simulation techniques in close collaboration with experiments is inevitable in this field.
We review how first-principles methods can be combined in terms of a multicode approach to investigate challenging complex oxide structures with respect to their electronic and magnetic properties. The paper is organized as follows. We start with an overview of different first-principles methods relevant for our investigations, highlighting the main concepts of each method together with its strengths and weaknesses. A sufficient treatment of the strong correlations appearing in oxides is necessary for a precise calculation of physical properties of oxides. Two schemes, the self-interaction corrections (SICs) and Hubbard corrections, will be discussed in more detail. Furthermore, we review how disorder can be described, which plays an important role in oxides. It can appear as disorder in a sublattice of a solid solution, as chemical disorder at interfaces, or in terms of intrinsic defects and impurities. Here, we introduce and compare two different approaches, the coherent potential approximation (CPA) to be used in connection to Green's function methods and the special quasirandom structures (SQS), which can be used in connection with efficient basis set methods allowing us to handle large supercells.
In addition to electronic properties we will explain the determination of magnetic ground state structures by first-principles methods together with the transition between magnetic structures in dependence on temperature. While the magnetic ground state structure can be found by comparison of the total energies of different magnetic structures, the temperature dependence will be investigated by means of an effective classical Heisenberg model. We review first-principles methods to determine Heisenberg exchange parameters as well as the implementation of the Monte Carlo method used to study the magnetic ground state and temperature-dependent behavior. In the last section, we review several applications and combinations of the methods on the example of relevant functional oxides. A summary concludes the paper. electronic structure requires nevertheless a nontrivial treatment of a many-body problem. The corresponding many-body wave function of a quantum mechanical system is a function of the position and spin of all particles, which results in a huge number of degrees of freedom. The calculation of such a wave function is possible only for very small systems consisting of a few atoms. Therefore, a quantum mechanical description of many-particle systems necessarily requires suitable approximations.

Basics of Density Functional Theory
A prominent simplification is given by the framework of the DFT, a first-principles approach that has been applied successfully for more than 50 years in simulations throughout solid state physics, nuclear physics, and chemistry. The DFT is based on the two Hohenberg and Kohn theorems, which state that all ground state properties of a quantum mechanical system can be uniquely described via its particle density and that a related and unique ground state exists. [15] According to the Hohenberg-Kohn theorems, the total energy of a system is a unique functional of an external potential vðrÞ caused by the underlying lattice of atomic nuclei and sources other than the solid itself which is determined, except for a trivial additive constant, by the electron density ρðrÞ. The electronic density ρðrÞ determines uniquely the ground state and all other electronic properties of the system, where T½ρ and V ee ½ρ are functionals of the kinetic energy and the electron-electron interaction, respectively. In Equation (1), the total energy is exact for any quantum mechanical system, but a closed expression for T½ρ and V ee ½ρ as functionals of ρðrÞ is unknown and, therefore, the equation cannot be directly exploited. Kohn and Sham suggested to represent the electronic density via electronic wave functions in the independent particle approximation. [16] In this approach, the total energy functional is given by (this and the following equations are expressed in Rydberg units) where T s ½ρ is the single-particle kinetic energy functional of the independent electrons, the second term is the same as in Equation (1), the third term is the Hartree or Coulomb energy functional, and the last term is the so-called exchange-correlation energy functional E xc ½ρ. The latter contains the difference between T½ρ and T s ½ρ and those contributions of V ee ½ρ which go beyond the electrostatic interaction of charge densities. This term is presumably small and must be approximated. The electronic density in the independent particle approximation is given by where the ψ i ðrÞ represent the N lowest occupied eigenstates of the one-electron Kohn-Sham Hamiltonian and are solutions of the corresponding Kohn-Sham equation ðÀ∇ 2 þ v eff ðrÞ À ε i Þψ i ðrÞ ¼ 0 (4) with an effective Kohn-Sham potential v eff ðrÞ ¼ vðrÞ þ 2 Z dr 0 ρðr 0 Þ jr 0 À rj |fflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflffl ffl} v H ðrÞ þ v xc ðrÞ (5) Here, the first term is the external potential, the second term is the Coulomb or Hartree potential v H ðrÞ, and the last term is the exchange-correlation potential defined as v xc ðrÞ ¼ δE xc ½ρ δρðrÞ (6) which must be approximated together with E xc ½ρ. The first efficient approximation was introduced by Kohn and Sham, who restrict E xc ½ρ only to the local density at r and approximate it with the uniform electron gas formula, [16] known as the local density approximation (LDA) or the local spin density approximation (LSDA) for magnetic systems (only the term LDA is used in the following independently of the particular system). Nowadays there exist various efficient parameterizations of E xc ½ρ. [17,18] Throughout our examples and studies discussed subsequently, we always chose the exchange-correlation functional best suited to the posed problem and most reasonable in terms of computing times. Where the parameterization of Perdew and Wang [19] for the LDA failed to describe the materials properties well enough, we adapted mostly the generalized gradient approximation (GGA) from Perdew-Burke-Ernzerhof (PBE). [20] In other cases, the especially for solids revised PBE functional PBEsol was applied [21] -we point out explicitly, when other functionals were used. In case of magnetic systems, the respective spin-resolved functionals were applied.
Altogether, Equation (4) reduces the complex many-body problem to that of one electron interacting with the effective field of all other electrons and the nuclei, expressed by v eff ðrÞ in Equation (5). The latter, however, also depends on the electronic charge density and the final solution is found by a self-consistent iteration of the so-called Kohn-Sham Equations (3)-(5). The total energy is calculated from Equation (2). Based on this iteration procedure, various numerical methods have been developed over the past decades being successful in describing ground state properties of materials systems. They can be divided into a few basic categories depending on the approach to solve the Kohn-Sham equations. In the following, we will focus on two methods: first, those based on basis sets and, second, the multiple scattering theory applied to the electronic structure of solid materials, which was pioneered by Korringa, Kohn and Rostoker (KKR) [22,23] and developed later to the Green's function method. [24,25] Both names, KKR and Green's function method, are often used synonymously. [25]

Basis Set Methods
The solution of the Kohn-Sham-Schrödinger Equation (4) can be obtained by expanding the Kohn-Sham wave functions in terms of a basis set. As a result, the differential equation is transformed into a system of linear algebraic equations for the expansion coefficients, which can be solved numerically. Although the general formalism remains the same, additional conditions and approximations might emerge from the choice of the respective basis. A prominent example is the introduction of so-called pseudopotentials that allow us to use efficiently a set of plane waves instead of strongly oscillating atomic wave functions. Such an approach is implemented in several program packages, such as VASP, [26,27] ABINIT, [28] and Quantum ESPRESSO. [29] In addition to plane waves, the linear combination of atomic orbitals (LCAO) offers another well-established approach, which is used in program codes such as CRYSTAL, [30] FPLO, [31] and SIESTA. [32] Furthermore, both approaches can be combined, e.g., in the full-potential augmented plane wave methods combining atomic orbitals and plane waves: WIEN2K, [33] FLEUR, [34] and ELK. [35] Another combination of atomic orbitals with partial waves leads to the linear muffin-tin orbital method (LMTO) [36,37] or the augmented spherical wave method. [38,39] In the following, we give an overview of the basic concepts of pseudopotential methods and LMTO methods.

Pseudopotentials
The modern pseudopotential and projector augmented plane wave (PAW) methods have the orthogonal plane wave method (OPW) [40] as a predecessor. The basic idea is to transform the original problem to a form that operates with smoothly varying functions, i.e., removing strongly oscillating parts of the wave function close to the nuclei. To briefly outline the main concept, we recall the very early ideas by Phillips and Kleinman. [41] Therein, the exact wave function of a valence electron ψ ν can be expressed as  (4)). The energies E c are the energy levels of the core states and the potential part in H PS is given by The nonlocal operator V R represents a repulsive potential and cancels the attractive potential from the cores, resulting in a smoothly varying pseudopotential V PS . In addition, the states of the core electrons are treated as constant-the so-called frozen core approximation.
Over the years several important schemes to construct pseudopotentials have been developed, such as norm-conserving pseudopotentials, [42] ultrasoft pseudopotentials, [43] and the PAW method, [44] to name only the most important. These approaches are implemented in packages such as VASP, ABINIT, or Quantum ESPRESSO.
The use of pseudopotentials is the basis of an efficient representation of the crystal wave function. A big advantage of the here-discussed codes is that they allow a detailed investigation of the real structure of complex bulk materials, defects, surfaces, or interfaces. Instead of using experimental information, which is not always available or precise, such methods allow for determining the ground state crystal structure in a theoretical consistent way by structural optimization.
To find the correct structure representing the energy minimum from an initial guess in agreement with all symmetry requirements, one has to minimize the forces acting on the atoms. The calculation of forces is based on the Hellmann-Feynman theorem (force theorem), where the force acting on an atom at position R is given as a derivative of the total energy E by [45] with the external potential vðrÞ in Equation (5) and the ionic energy E I . The latter can be efficiently calculated using a DFT approach with an appropriate exchange-correlation potential. However, the evaluation of forces in Equation (10) is an intricate task, which cannot be implemented efficiently in all DFT methods. For example, the Green's function method or the LMTO-ASA method described later cannot provide accurate calculations of forces because of their potential representation and the partial wave expansion of the wave function. Therefore, plane wave methods are often used to generate ground state structural information, which serves as input for other methods, as presented for our multicode approach later, where VASP is one of the basic tools.

Linear Muffin-Tin Orbitals
The LMTO method can be derived from partial wave methods by energy linearization. Generally, the atomic sphere approximation (ASA) is used, i.e., the potential is represented by slightly overlapping, space-filling spheres around the nuclei (LMTO-ASA), although there exist several efficient full-potential implementations of the LMTO method. [46][47][48] In LMTO-ASA, the potential inside the sphere is spherical symmetric. The solution of the single-site problem in the sphere at R with angular momentum L ¼ ðl, mÞ and energy E is given by φ R,L ðr, EÞ ¼ φ R,l ðr, EÞY L ðrÞ. The total wave function ψðrÞ is built from so-called muffin-tin orbitals defined inside and outside the sphere as www.advancedsciencenews.com www.pss-b.com The functions N Rl ðEÞ and P Rl ðEÞ are normalization and potential functions, while J L ðrÞ and K L ðrÞ are solutions of the Laplace equation [J L ðrÞ ¼ J l ðrÞY L ðrÞ, K L ðrÞ ¼ K l ðrÞY L ðrÞ]. This formulation leads to a KKR-like secular equation that separates the properties of the individual atoms from the structure of the system, but has a nonlinear energy dependence. The key step is to linearize the energy dependence of the radial function by a Taylor series expansion around some fixed energy E 0 ¼ E ν,Rl , which represents usually the center of the occupied part of the R and l projected density of states of the valence electrons Here,φ R,l ðrÞ is the energy derivative of φ R,l ðr, EÞ at E 0 . This linearization restores an eigenvalue problem. Lattice Fourier transformations are used to take translational symmetry into account. Despite the limitations due to the potential approximation, the LMTO method allows for the construction of Green's functions, the treatment of disorder, or the extension into a relativistic theory. [49] In addition, next generation muffin-tin methods are also applied successfully to various materials research problems, e.g., the exact muffin-tin orbital method (EMTO). [50] In our studies of oxide systems, we used the fully relativistic LMTO band structure method for simulations of magneto-optics, X-ray absorption spectroscopy (XAS), and X-ray magnetic circular dichroism (XMCD) data (see Section 6.3). It applies an implementation, which uses four-component basis functions constructed by solving the Dirac equation inside an atomic sphere. [51,52] Further details of the method can be found in previous studies. [49,53,54]

Green's Function Method
In contrast, the solution of the Kohn-Sham Equation (4) can be approached also in terms of the multiple scattering theory, which leads in a natural way to the Green's function method. [24,25,[55][56][57][58] This formulation provides flexibility and wide applicability, because the Green's function itself contains all information of the system and can be straightforwardly used to calculate the spectral properties of the Kohn-Sham Hamiltonian, the charge density of the system, or any other observable. [22,23] Without providing a complete list, there are various successful realizations of the Green's function approach available, such as the Munich SPRKKR program package, [58] the Jülich massively parallel Green's function method for large-scale systems, [59] the LSMS package, [60] AkaiKKR, [61] the disordered local moment method (DLM), [62] or HUTSEPOT. [63,64] The latter is the method of choice for our multicode approach, as it is designed for a multitude of problems in materials research from bulk materials, surfaces, interfaces, or real space clusters. The code development started several years ago with the calculation of semirelativistic angle-resolved photoemission spectra. [65] In general, the one-electron Kohn-Sham Green's function of a complex energy E can be constructed in real space in terms of a complete set of orthonormal wave functions, ψ i ðrÞ, the eigenfunctions of a Kohn-Sham Hamiltonian, and corresponding eigenenergies ε i Gðr, r 0 ; EÞ ¼ Alternatively, the Green's function Equation (15) can be found as the resolvent of the Kohn-Sham Equation (4) by solving self-consistently ðÀ∇ 2 þ v eff ðrÞ À EÞGðr, r 0 ; EÞ ¼ Àδðr, r 0 Þ where v eff ðr) is the effective Kohn-Sham potential, as introduced in Equation (5). The corresponding Green's function can be constructed for any complex energy E. It determines the expectation value of an observableÔ by where f ðEÞ is the Fermi-Dirac distribution function. In particular, the electronic density in Equation (3) can be also calculated using the Green's function as Within the multiple scattering theory, a crystal under consideration has to be decomposed into distinct atomic regions. Therefore, we use in HUTSEPOT various shape approximations for the effective Kohn-Sham potential in Equation (5). Assuming a nonoverlapping spherical potential at each atomic site with a constant interstitial region utilizes the so-called muffin-tin approximation (MTA). This method is simple to implement and is very efficient for metallic close-packed systems. Another spherical approach is given by the ASA, where the volume of all atomic spheres in the unit cell is chosen to be equal to the volume of the Wigner-Seitz cell. This method provides a better description of open systems than the MTA, although it suffers from an artificial scattering of electrons because of overlapping atomic spheres. The ASA Green's function can be improved by inserting so-called empty spheres (with a nuclear charge Z ¼ 0) into the interstitial region. However, the error introduced by the overlap of the atomic potentials remains, although it can be reduced efficiently using a proper choice of atomic and empty spheres. The accuracy of total energy calculations can be substantially improved using a full-charge density representation in the integrals entering Equation (2). [66] In this method the charge density (Equation (18)) is expanded into spherical harmonics Y L ðrÞ ¼ Y L ðθ, ϕÞ inside the Wigner-Seitz cells as where L ≡ ðl, mÞ. The integrals in Equation (2) are calculated using a shape function representing the underlying Wigner-Seitz cell. In contrast, the real shape of the Wigner-Seitz cell can be taken into account by a full-potential approach, [57,61,67,68] in which the effective potential is decomposed with respect to spherical harmonics v eff ðrÞ ¼ Based on this full-potential approximation and the real-space multiple scattering theory, the Green's function is expressed for any arrangement of atoms in terms of the scattering path operator τ nm ðEÞ Gðr n , r 0 m ; EÞ ¼ The building blocks of Green's function are the regular solutions, Z L ðr; EÞ, and the irregular solutions, J L ðr; EÞ, of the radial Schrödinger equation at the given (complex) energy E (e.g., Kohn-Sham Equation (4)), which are matched to spherical Bessel and Hankel functions outside the potential range (r ≥ S). [69] The scattering path operator describes all scattering events between different lattice sites and is given in the angular momentum representation. For general electronic systems, the τ-matrix is implicitly given in terms of the single-site scattering t-matrices of the different atoms and the structure constants gðEÞ, representing the free-electron Green's function, and can be found from the matrix equation The t-matrix can be determined from the normalization conditions of the wave functions and describes the scattering of an electron at the atomic potential. [69] Equation (22) is the main quantity of the Green's function method. It yields a complete separation of the atomic properties at different sites of a material, expressed in the scattering matrices, from the structural aspects, embodied in the structure constants gðEÞ of the underlying lattice. The structure constants gðEÞ can be obtained in real space and generalized to any symmetric lattice by a corresponding Fourier transformation according to the symmetry of the problem.
Moreover, the Green's function formalism can be naturally extended to incorporate also relativistic effects by solving a singlesite Kohn-Sham Dirac equation instead of the Kohn-Sham Schrödinger Equation (4). [68] This approach was implemented also in HUTSEPOT within the full-potential approximation. [70]

Group Theory and GTPack
Symmetry and symmetry breaking, as general concepts of nature, play a significant role also in oxide materials. While the crystalline symmetry of oxides determines the degeneracy of electronic, magnonic, or phononic bands as well as the degeneracy and magnetization of local defect states, symmetry breaking corresponds to the underlying mechanism behind, e.g., interaction instabilities driving the system into various ordered phases of matter such as magnetic, orbital, or charge ordering as well as superconductivity. To provide a powerful tool to complement first-principles investigations with symmetry analyses, we have implemented the Mathematica group theory package GTPack. [71] GTPack contains more than 200 additional modules to the Mathematica core language, covering several fields of application such as basic group theory and tables of point and space groups, representation theory, crystal field theory, electronic structure calculations in the tight-binding and pseudopotential approximations, photonic band structure calculations, and Landau-Ginzburg-Wilson theory of phase transitions. Some functions allow an import of data from firstprinciples calculations for further investigations. The construction of tight-binding models helps to understand better the complex physics, using input parameters from the firstprinciples calculations.
Among others, GTPack was successfully applied for the classification of superconducting gap functions including odd-frequency superconducting states as they might arise in strongly correlated layered cuprates [72] or the analysis of topologically protected Hund-nodal-line semimetals, [73] which can be realized in double-perovskite manganites in the twisted magnetic phase. Details about the implementation of GTPack and a comprehensive introduction to group theory are given by Geilhufe and Hergert. [74,75]

Correlation Corrections
As many oxides contain transition metal or rare earth elements, strongly localized d-and f-electrons arise, which exhibit significant Coulomb interactions compared to the order of magnitude of their kinetic energy. These localized states in the so-called strongly correlated materials represent a challenge for firstprinciples investigations, because conventional LDA and GGA exchange-correlation functionals typically fail to describe correctly their electronic structure. [12,13] To overcome this limitation, several correlation corrections have been developed over the past decades (see previous studies [76,77] for a detailed overview).
The most popular are the SIC method, [78] the LDAþU approach, [79][80][81] and the extension of the last-the dynamical mean field theory (DMFT). [82,83] While the first two can well reproduce ground state properties of strongly correlated materials, the DMFT can describe adequately excited state properties and other many-body effects, e.g., metal-to-insulator transitions in TMOs. [84][85][86][87] The DMFT uses the spatial localization of the Coulomb interaction to map the many-body lattice problem to a local but time-or energy-dependent mean field, [83] hence the name dynamical mean-field theory, where at its basis a singleimpurity Anderson model has to be solved self-consistently. [76] In many cases such an impurity model gives accurate results but is usually solved with computationally expensive quantum Monte Carlo methods. Therefore, the SIC and LDAþU approaches have remained dominant in many theoretical studies and they are most crucial for our investigations on complex oxides as well (see Section 6). In particular, we implemented a numerically simpler scheme of SIC in terms of a local SIC method in HUTSEPOT, [88] which we outline in Section 3.2. But before, we briefly review the necessity of the SIC and underlying ideas.

Self-Interaction Correction
The need for a SIC follows directly from the DFT scheme and the approximations to the exchange-correlation functional, because the latter is not strictly self-interaction-free, whereas the exact density functional is, as observed by Perdew and Zunger [78] as well as other authors before (see, e.g., previous studies [89,90] ). This self-interaction term should vanish exactly, as it is realized in Hartree-Fock theory, for a single, fully occupied orbital, expressed here in terms of a density ρ ασ with quantum number α and spin σ [78] where E H is the Hartree energy v H ðrÞ realized from the second contribution in Equation (5) and E exact xc is the unknown exact exchange-correlation energy in Equation (6).
However, the cancellation is incomplete and Perdew and Zunger [78] proposed an orbital-dependent self-interactioncorrected energy functional which is in principle valid for all approximate E xc which do not satisfy Equation (23). The tilde forẼ approx denotes that this energy functional has basically the same form as its equivalent in the Kohn-Sham theory but the orbitals in Equation (24) do instead minimize the self-interaction corrected energy functional. [91] In particular for the LSDA, Equation (24) leads to an orbital-dependent effective SIC potential seen by an electron in orbital ψ ασ v SICÀLSDA eff ,ασ with the potentials in Equation (5).
One possible implementation of a generalized version of the scheme presented by Perdew and Zunger [78] (referred to as full SIC subsequently) is realized in the LMTO-ASA band-structure method, [36,92,93] where the complete implementation of the SIC-LSDA formalism is provided in previous studies. [94,95] One advantage of this formalism is the possibility to study the valency of ions. The total energy minimization is not unique anymore but depends on the localized orbitals as well. When the lowest configuration with N SIC localized states is found, the resulting valency is with the atomic number Z and the number of core states N core . The valency is successfully reproduced for various materials from rare-earth elements, actinides, over TMOs, to materials exhibiting the colossal magnetoresistance effect. [96][97][98][99][100] In contrast, the full SIC-LSDA scheme is rather time consuming because it involves a band-structure problem. The scheme makes it necessary to repeat transformations from reciprocal space to real space, transform Bloch functions to Wannier functions, and evaluate the SIC potential followed by a back transformation to reciprocal space. [88,91]

Local Self-Interaction Corrections
To minimize computational effort and take advantage of the multiple scattering formalism, we introduced and implemented the so-called local SIC (LSIC) to HUTSEPOT. [88,91] It utilizes the very sharp resonances of the scattering phase shifts of bound states, associated with a large Wigner-delay time at the considered site. Localized valence states show an abrupt jump of π in the generalized complex phase shift similar to the core electron states but at positive energies. [91] This complex phase δ i l ðϵÞ contributes to Equation (22) as the t-matrix with κ being ffiffi ffi ϵ p , the site index being i, and the angular momentum quantum number being l.
As before, the charge density is needed for the self-consistent cycle and we can in accordance with Equation (18) define an orbital-resolved corrected density whereas the angular momentum contributionsL ≡ ðl,mÞ result from an unitary transformation in order to adapt the spherical symmetry. [91] The integration in Equation (28) runs over an energy contour encompassing the whole valence band. [91] The resulting density is used for self-interaction-corrected potentials expressed as The latter potential defines the single-scattering matrices where some orbitals marked withL c σ c are corrected and the t-matrix element t iL σ is obtained from the uncorrected LSDA potential. [91] Finally, these t-matrices are put into Equation (22) to obtain the self-interaction-corrected scattering-path operator, which is then used in the self-consistent cycle again to calculate the total energy (Equation (2)). The latter is also orbital dependent and allows again the study of different valence states via Equation (26).
The LSIC scheme was demonstrated and compared with the full SIC approach implemented in the LMTO code at first for NiO as a prototypical strong correlated material. [91] NiO is one compound of the 3d TMO series (MnO, FeO, CoO, NiO, and CuO) and appears as a comparative example at several places in this work. Hence, we briefly recollect the basic characteristics of the TMO. All of them crystallize in the rock-salt lattice structure (see Section 4 and 5 for a structural sketch). They are wide-gap Mott-Hubbard or charge transfer antiferromagnetic (AFM) insulators and order AFM as type 2 (AFII): Two planes in [111] direction always have opposite spin and create two sublattices. The TMOs are one example where the convential LSDA und GGA exchange-correlation functionals fail and correlation corrections improve the electronic structure, resulting in magnetic moments and bandgaps in agreement with experiments. [80,[84][85][86][87][91][92][93][102][103][104][105] In particular for NiO, the LSIC recovers well the important features of the density of states (DOS) calculated also with the full SIC approach (Figure 1a). The corrected Ni d states are shifted correctly to the bottom of the valence band and the Ni3d-O 2p hybridization is reduced. The comparison to experimental valence band photoemission spectra reveals also a better agreement than the DOS calculated with the conventional LSDA ( Figure 1). Even though the Ni d states are slightly higher in energy than in the full SIC result, the bandgap and the O 2p states close to the Fermi energy are almost similar. The provided bandgap of 3.7 eV does not exactly match with the experimentally measured value of 4.3 eV but is a large improvement with respect to the DOS obtained with LSDA, which is even metallic with the Ni 3d states localized in the upper part of the valence band ( Figure 1b). We note that the Fermi energy obtained from the numerical calculation lies inside the bandgap but cannot be directly related to the experimental spectra. Although in calculations the Fermi level can be placed arbitrarily between the valence and conductance bands, the position of the Fermi level in experiments can be affected by any kind of doping and other imperfections. Hence, we shifted in Figure 1-3 the calculated DOS by 1 eV with respect to the original publications.
For the LSIC results, we have, in contrast to conventional selfconsistent calculations, at first determined the ground state configuration and valency, because the total energy and the potential in Equation (29) are orbital dependent. These are identified by the lowest total energy in dependence of different selfinteraction-corrected orbitals-here, the 3d states of Ni. It turned out that the resulting 2 þ state, with electrons in the five majority Ni 3d and three minority Ni t 2g states, fulfills indeed the first Hund's rule (maximizing the spin moment). [91] The corrected states become very sharp in the DOS and are shifted by % 10 eV down in energy, whereas the uncorrected states are shifted to higher energies ( Figure 2a). The well-localized character of the considered Ni 3d states can be also verified via the phase shifts, as discussed previously. Whereas the corrected states have a very sharp resonance, the uncorrected states or the states considered in the LSDA calculation vary continuously ( Figure 2b).

Correlation Correction via Hubbard U Parameter
Another approach to correct the description of localized states is based on the assumption that the strongly localized states can be described via an on-site atomic interaction in terms of the Hubbard Hamiltonian. [106] Anisimov and coworkers [79,107,108] suggested to include this interaction at the DFT level via the Hubbard U which can be understood as the Coulomb energy required to place two electrons at the same site, i.e., for d or f states. In those works, the authors introduced the LDAþU method, but the same concept (a) (b) Figure 1. The NiO total density of states calculated with a) the full SIC approach implemented in the LMTO-ASA [94,95] (blue line) and LSIC (black line) and b) conventional LSDA (gray line). The latter two were obtained with HUTSEPOT. The right axis in (b) is related to the intensities of simulated and experimental valence band photoemission spectra for different photon energies from Taguchi et al. [106] Adapted with permission. [91] Copyright 2009, Publisher IOP Publishing Ltd. can be applied in principle to a generic approximate correlationexchange functional. [81] Hence, the term 'GGAþU' also appears subsequently with respect to our examples in Section 6. The LDAþU corrects the original DFT total energy of the localized states based on a Hubbard-like term E U , while the other states are treated normally [81] -a similar idea as SIC in Equation (24). However, an additional double-counting term E dc has to be subtracted to eliminate the part of interaction energy which is already included in E U . [81] Although nowadays several different ways of Hubbard U corrections are implemented in many first-principles methods, we will use the most basic formulation to outline the underlying concept. [79,108,109] Let us assume that all d electrons of a specific ion have roughly the same kinetic energy T d and feel the same Coulomb repulsion U. The total energy of a shell with n electrons is E n ¼ nT d þ nðn À 1ÞU=2 with the excitation spectrum tions as a function of the number of d electrons N is in such a case E ¼ NðN À 1ÞU=2. [108] In fact, the latter represents the double counting term E dc in the most general form. It has to be subtracted from the total energy (Equation (2)) when E U is added [108] E with the orbital occupation numbers n α . The orbital energies are derivatives of Equation (32) and the effective potential in the Kohn-Sham Equation (4) becomes orbital dependent Both of the previous equations reveal the main features of the LDAþU approach. The total energy of occupied orbitals (n α ¼ 1) is according to Equation (33) shifted down by ÀU=2, whereas the energy of unoccupied orbitals (n i ¼ 0) is shifted up by þU=2. This opens a bandgap in case of the Mott insulators. At the same time, the correction to the potential restores the discontinuous behavior of the one-electron potential of the exact DFT. [81,108] In modern implementations of the LDAþU, the exchange interaction and a nonsphericity of the on-site Coulomb interaction are taken into account by a rotationally invariant formulation [110] E LDAþU ½ρðrÞ, fρ σ g ¼ E LDA ½ρ σ ðrÞþ þ E U ðfρ σ gÞ À E dc ðfρ σ gÞ (35) where E U ðfρ σ gÞ and E dc ðfρ σ gÞ depend now on elements fρg of the density matrixρ σ mm 0 (σspin direction, mmagnetic quantum number). The latter consists in the general case of spin-diagonal and spin-off-diagonal terms and can be calculated, e.g., from Green's function.
A simplified rotational invariant formulation of the LDAþU is implemented in both program codes, VASP and HUTSEPOT, and given explicitly as [80] where U and J are spherically averaged matrix elements of the screened Coulomb and the exchange interactions, respectively. Only the difference U eff ¼ U À J as an effective U parameter is relevant in Equation (36). The latter was explicitly used to obtain agreement between measured and calculated electron energy loss spectra (EELS) for NiO with U eff ¼ 5.3 eV. [80] Other calculations, with LSDA alone or with larger U eff parameters, did not agree with the experimental results. Therefore, we used a similar U eff to compare the DOS of NiO with the ones calculated with LSIC and LSDA ( Figure 3). The LDAþU recovers as well as the LSIC the insulating character of the Mott insulator NiO, which is not available within LSDA. However, the bandgap of 2.7 eV obtained from LDAþU is smaller than the one obtained from LSIC. The Ni 3d states are also shifted down in energy, forming broader valence bands due to a hybridization of Ni 3d and O p states but they are not as separated as it could be expected from the model in Taguchi et al. [106] (see Figure 1). In summary, the LDAþU reproduces for NiO to some extent properties of excited states, e.g., EELS, [80] whereas the LSIC offers a parameter-free method to describe its ground state properties. Nevertheless, the most accurate metal-to-insulator transition temperature is obtained with DMFT. [87]

Advantages and Disadvantages of the Two Methods
The SIC and LDAþU approaches improve the description of structural, electronic, and magnetic properties of strongly correlated systems in comparison with conventional LDA or GGA functionals. However, it is not directly possible to compare the two. They are based on different ideas to correct the total energy in Equation (1): SIC acts at the level of the exchangecorrelation functional and subtracts the Hartree and LDA total Figure 3. Density of states for NiO in the AFII structure (gray curve) calculated using the LDA þ U method within HUTSEPOT (U eff ¼ 5.3 eV motivated from Dudarev et al. [80] ). The structure was taken from experimental measurements. The estimated bandgap value is 2.7 eV. For comparison, the simulated and experimental valence band photoemission spectra for different photon energies from Taguchi et al. [106] and results for LSIC are replotted from Figure 1a. www.advancedsciencenews.com www.pss-b.com energy for a specific charge density from the whole density, whereas LDAþU introduces a physically motivated parameter. Nevertheless, both methods are highly needed for the study of oxides, because those include quite often strongly localized d or f electrons (as demonstrated with the examples in Section 6). The numerical implementation of the SIC scheme is, most importantly, a first-principles approach, as the implementation of SIC does not involve any additional physical parameters. The determination of the ion valency becomes very crucial for magnetic oxides with a low, high, or intermediate spin state as, e.g., in SrCoO 3 , [110] or oxides with defects like oxygen vacancies. At the same time, the SIC method is designed mainly for the description of ground state properties and cannot describe excited state properties very well. Here, the combination of SIC with Slater transition approaches [111] offers an approximation of these characteristics. [91] However, systems with weaker correlation effects might suffer from an overlocalization of the orbitals within SIC (e.g., in Sr 2 FeMoO 6 , Section 6.2), where the LDAþU can perform better. By comparing, numerical results for different U parameters with experimental results for, e.g., photoemission spectra, XAS, and XMCD, Curie temperature, or spin waves, one can estimate the impact of electronic correlation to those quantities or find a good approximation for the underlying electronic structure of the considered material. In addition to this semi-empirical approach, there exist several methods to calculate the Hubbard U parameter from first-principles, e.g., Slater's transition state technique, linear response, [112] from constrained random phase approximation calculations. [113] However, for the materials considered in our studies either no U parameters were calculated or they were strongly overestimated. One example for an overestimated U parameter is found for SrCoO 3 . [114,115] The calculated U ¼ 10.83 eV and J ¼ 0.76 eV parameters are far too large to describe well the magnetic transition temperature of % 300 K.
Other advantages are implicitly included. Whereas the implementation of LSIC into HUTSEPOT offers the features of the multiple scattering theory, as disordered systems treated with CPA or NLCPA (see Section 4), the simple formulation of the LDAþU approach allows the calculation of energy derivatives. [81] Forces as in Equation (10), stresses, or dynamical matrices are most crucial to find the ground state structure of the considered oxides.

Treatment of Disorder
Almost all numerical first-principles methods, as discussed in Section 2, are designed to treat 3D periodic systems. This allows us to exploit Bloch's theorem to reduce the numerical effort. Low-dimensional systems, such as surfaces and clusters, can be treated in a supercell approach, to restore the translational symmetry at the expense of increased numerical burden. Realspace formulations do not rely on Bloch's theorem and are of advantage if systems such as quasicrystals are studied. [116][117][118] Due to chemical disorder in solid solutions, e.g., A x B 1Àx , the translational symmetry is broken. Although the concentration x of the species is known, the microscopic distribution of the atoms in the sample remains hidden. Two principal approaches to perform a configuration average are possible: the CPA in the framework of multiple scattering theory constructs an effective medium, having the advantage to be lattice periodic again [119][120][121] ; the construction of SQS mimics the correlation functions of an infinite substitutional alloy in a certain supercell. [122,123] Subsequently, we review the general concepts and implementations of both methods and give respective examples to discuss and compare their performance.

Coherent Potential Approximation
The general idea of the CPA is to construct an auxiliary medium, characterized by a single-site t-matrix t C which takes into account implicitly all configurations of the disordered alloy A x B 1Àx . The construction of this effective medium in multiple scattering theory is based on the single-site t-matrices t A ðEÞ and t B ðEÞ, describing the scattering properties of the individual atoms A, B and the scattering path operator, as defined in Equation (22).
For the site-diagonal quantity τ C , we get for the binary alloy i.e., τ C is found from the concentration-weighted scattering path operators τ A , τ B describing the scattering properties of A,B embedded in the effective medium (see Figure 4). As an example τ A can be obtained from The scattering path operator of the coherent medium can be evaluated for 3D periodic solids from where Ω BZ is the volume of the Brillouin Zone (BZ) and gðk; EÞ represents the structural Green's function (see Equation (22)) in momentum space representation. The relation (39) between t C and τ C together with Equation (37) and (38) allows for a selfconsistent determination of the scattering properties of the effective medium. The scheme can be easily generalized to N components each with concentration c i and i ¼ 1,⋯,N. A new approximation t new C can be generated from the previous step by By means of Equation (40), the new scattering path operator of the effective medium is calculated. The CPA condition (Equation (37)) is an additional self-consistency requirement on top of charge or potential self-consistency.
The CPA equations can be easily generalized for any boundary conditions. In contrast to many simple single-site theories such The sites labeled C are characterized by a coherent potential, while sites labeled A and B represent the impurity potentials. The impurity concentrations at both sites are not independent: www.advancedsciencenews.com www.pss-b.com as virtual crystal or average t-matrix approximations, the CPA is exact in both, the weak-scattering and the narrow-band, limits. It has been successfully applied to basic problems in statistical physics, [124] electronic structure calculations of bulk and lowdimensional materials, [125][126][127] and so on. However, the CPA remains a single-site approximation and environmental effects on scattering properties are neglected, except on average. This limitation does not allow us to investigate fluctuations around the CPA average and to elucidate the influence of atomic short-range order. Such multisite effects can be systematically taken into account using a nonlocal coherent potential approximation (NLCPA). This theory was recently derived within the KKR framework. [128][129][130][131][132] The books of Gonis [25,56] and Turek et al. [49] present more basic details and implementation-related discussion of the treatment of disorder.

Special Quasirandom Structures
In a solid solution A x B 1Àx , the local structural relaxations around an individual A or B atom play an important role, even if the experimental techniques probe local properties only on average. An obvious drawback of the CPA is that local relaxations due to the microscopic distribution of the constituents cannot be taken into account. In contrast, a binary system with N sites leads to n c ¼ 2 N different configurations σ, which have to be relaxed and averaged; i.e., to get the average value P of the physical property P one has to evaluate where wðσÞ denotes the probability of a certain configuration in the ensemble. The construction of SQS is an approximation, pioneered by Zunger et al., [122,123] to handle the large configuration space. It is based on the cluster expansion formalism. [133][134][135] The method is based on the fact that contributions from distant neighbor interactions to the total energy become negligible beyond a certain range. A configuration is discretized in figures f ¼ ðk, mÞ consisting of k vertices (pair, triplet, tetrahedron) and separated by an mth neighbor distance. This decomposition is used to mimic the correlation functions of the solid solution as closely as possible. Details of the method can be found in previous studies. [122,123,[136][137][138][139] The SQS method is designed for a rigorous consideration of binary, ternary, and multicomponent alloy systems by means of large supercells. An example for a disordered ternary system of MgO and ZnO is shown in Figure 5. A practical realization of the method is discussed by van de Walle et al. [138,139] The obvious advantage is that the normal DFT codes can be used. The relaxations in the supercell also include the microscopic structure around A and B atoms. The construction of the SQS is separated from the electronic structure calculation itself.

Comparison of the Two Methods
The topic of disordered systems is still an interesting and heavily investigated problem in materials sciences. CPA and SQS compare very well for alloys, as demonstrated in several works for Fe-Cr, Ti-Nb, or Ti-V, [141][142][143] and Ti-Al or Cu-Au alloys. [144,145] Both methods are applied to different kinds of systems, including oxides [146][147][148] and high-entropy alloys. [149][150][151][152] Therefore, our studies usually apply both methods. For example, Maznichenko et al. [125] investigated structural phase transitions and fundamental bandgaps in Mg x Zn 1Àx O by means of CPA as implemented in HUTSEPOT. According to Figure 6 three structural phases appear in the phase diagram, namely, rock-salt (RS), hexagonal (HX), and wurtzite (WZ) structure. In addition to the treatment of disorder, correlation corrections to the Zn 3d electrons in terms of the SIC scheme were included.
In the following, both methods to treat disorder, CPA and SQS, are compared on the example of Mg x Zn 1Àx O. A concentration of x ¼ 0.75, i.e., a solid solution in the RS phase with the corresponding equilibrium volume, was selected for comparison. The SQS, constructed for this concentration, is shown in Figure 5. The SQS was constructed by means of the Alloy Theoretical Automated Toolkit (ATAT). [138] The supercell was built with 64 atoms starting in RS structure, which contains only two atoms. To sufficiently mimic the correlation functions of a  For the actual calculations in KKR-CPA, as well as in VASP-SQS, the LDAþU scheme was used for correlation corrections. Figure 7 compares the DOS, where a good agreement of both calculations is obtained. It is quite common that a unique U eff value does not lead to identical results in different firstprinciples methods due to distinct approximations used in the codes; e.g., typically chosen values of U eff in a KKR setup are smaller than in a VASP setup. Because U eff is seen as an adjustable parameter, it can be used to achieve agreement between the methods in the basic electronic structure as a requirement to use the distinct strengths of the methods to calculate physical properties.

Magnetism
As many oxide materials are discussed as applications for magnetic devices, the accurate description of those magnetic properties is one of the most important aspects in oxide research. Many of them, such as spin densities, the magnetic ground state, and exchange interactions, are not easy or impossible to access by experimental means and need first-principles calculations for a better understanding. Perovskite manganites, e.g., show a series of interesting properties, including colossal magnetoresistance, coexistence of clusters with different phases, rich phase diagrams, and effects like unusual spin, charge, lattice, and orbital order. [153] The magnetic characteristics such as spin-polarized band structure, magnetic moments, spin densities, and magnetization can be obtained directly from the Kohn-Sham Equation (4) using an appropriate spin density functional. The ground state magnetic order can be determined from total energy calculations of various static spin configurations, although this method might be very demanding for systems with a nontrivial magnetic structure. It can be done more efficiently utilizing a noncollinear spin density functional, which enables a direct search of the ground state magnetic structure by simultaneous simulations of magnetization dynamics. [154,155] Temperature effects, on which magnetic properties are extremely sensitive, can be taken into account, e.g., by the disordered local moment picture, [62,156] which is implemented in HUTSEPOT. In this approach, a magnetic system is modeled as an array of microscopic magnetic moments of fixed magnitude but random orientation in terms of the CPA. This approach described, e.g., very successfully the variation of the magnetic transition temperature in SrCoO 3 with varying oxygen vacancy concentrations. [110] Other magnetic properties important for the study of oxides, e.g., magnetization dynamics and thermodynamics, including the determination of the spin-wave dispersion relation, or the transition temperature, can be studied with the Heisenberg model Hamiltonian Here, the unit vectors are e i ≡ S i =S i at site i, where S i is the localized spin moment with an always positive moment size S i . The J ij are the so-called exchange parameters. We adopt the convention that J ii ¼ 0. The exchange parameters represent one of our most important tools to understand the microscopic magnetic order and can either be determined from experiments or in particular calculated from first-principles methods.

Calculation of Exchange Interactions
In the literature, several approaches are presented to calculate the magnetic exchange interactions J ij . For example, the so-called frozen magnon approach is one possible technique allowing us to find directly the energies of different magnetic configurations. [154] In contrast, J ij can be also extracted from the knowledge of the static magnetic susceptibility. [157] In this work, however, we focus on two approaches being more compatible with our considered multiple scatting (Section 5.1.1) or pseudopotential methods (Section 5.1.2). We consider again the TMO series MnO, FeO, CoO, and NiO, as a simple oxide test case (see also Section 3), and describe the underlying ideas of the calculation of J ij by recapitulating their numerical results. [103] All of them are in particular good examples of the Heisenberg model (Equation (42)), as their exchange interactions are very short ranged and only two are relevant to capture their magnetic structure (see sketch in Figure 8).  The exchange constant J 1 represents the nearest neighbor interaction, while J 2 is a next-nearest superexchange interaction including the oxygen atom. [103]

Explicit Calculation from Scattering Quantities
By means of the so-called magnetic force theorem (MFT), one obtains an explicit equation to determine J ij . [158] The basic idea behind the MFT is to consider a small rotation of classical spins at two different lattice sites as perturbation and map the resulting energy change onto the Heisenberg model (Equation (42)). This energy change can be expressed via the Lloyd formula and thereby represented by using the Green's function. Hence, the J ij between two lattice sites depend only on the scattering path operator τ ij (already defined in Equation (22)) and on differences of the single site t-matrices Δ i ¼ t À1 i" À t À1 i# for the two spin directions at site i. Both quantities are easily accessed via the Green's function method and the J ij are calculated by Note, compared to the equation in Liechtenstein et al. [158] an average over two spin directions results in an additional factor of 1=2 in Equation (43).
In practice, one needs only the respective lattice structure and calculates the electronic structure via self-consistent iterations of the Kohn-Sham equations. Then the quantities of the multiple scattering theory determine Equation (43). Thus, it is easy to vary, e.g., the lattice constant a of the TMO and determine the dependence of the J ij (see Figure 9). The variation of lattice constants can be related to different pressures. [88] For example, the nearest neighbor exchange interaction J 1 is for NiO ferromagnetic (FM, positive) and almost constant at different pressures, whereas the other TMO show a strong increase of the absolute value of J 1 (Figure 9a). The compressed lattice causes a stronger overlap of the d orbitals with the p orbitals of oxygen. Thus, a stronger superexchange is obtained. A similar behavior is found for J 2 for all TMO (Figure 9b). Here, the NiO J 2 is strongly AFM, compensating the FM J 1 , while the AFM character reduces with the series CoO, FeO, and MnO. Other examples for calculated J ij with Equation (43) are given in Section 6.

Calculation from Total Energy Differences
The J ij can be also calculated from the total energy differences of several magnetically ordered structures. This is related to the idea of only perturbing the magnetic state by a small rotation and obtaining the infinitesimal total energy change as used for the MFT described previously. In general, any first-principles method can be used to derive the total energies, but we have a bias for the pseudopotential method VASP as being a reliable tool to determine DFT total energies. The arising energy differences are mapped again to the corresponding energy differences in the Heisenberg model, given by Equation (42). This results in a linear system of equations for the J ij and is a rather general approach. However, its feasibility depends on the range of the exchange parameters, which determines the number of magnetic structures that have to be considered.
For example for the TMO, only nearest and next-nearest neighbors J ij have to be taken into account, as shown in Figure 8. Thus, three magnetic structures are used: the FM and two AFM (AFI, AFII) configurations. The AFI structure is characterized by oppositely magnetized planes, which are stacked along the (001) direction, whereas those planes are stacked along the (111) direction in the AFII structure. [103] By counting the number of nearest and next-nearest neighbors for the TMO, the J ij are determined as The resulting J ij agree very well with those obtained from the MFT previously and experimental values ( Table 1). The complete discussion for all 3d-TMO is given by Fischer et al. [103] Another example of this energy difference method in oxides is provided by Ködderitzsch et al. [101] Equation (44) and (45) were applied to bulk NiO in the framework of the full SIC-LMTO. In addition, they were extended to study the NiO(100) surface.
(a) (b) Figure 9. The exchange parameters a) J 1 and b) J 2 calculated with HUTSEPOT by Equation (43) for different lattice constants a of TMO. The vertical lines mark the calculated equilibrium lattice constants. Adapted with permission. [88] Copyright 2009, American Physical Society. Table 1. Comparison of the exchange parameters (in meV) of NiO obtained with the energy difference method ΔE and MFT [103] with experimental results from Hutchings and Samuelsen. [159] The calculated magnetic moment for Ni was 1.68 μ B and the experimental values were adapted to match the same Heisenberg Hamiltonian (Equation (42)). [103] Exp.
MFT Due to the broken symmetry, more exchange coupling constants enter the model. [101] Recently, Ben Hamed et al. [148] investigated the magnetic phase diagram of Gd x Ca 1Àx MnO 3 using a model containing up to three exchange constants. A calculation of the exchange parameters by the MFT approach for the latter system reveals also that more distant interactions have to be taken into account. Hence, the model was extended to include eight exchange parameters calculated within the energy difference method with much more magnetic structures. [160] It results in an improvement of the calculated Néel temperature.

Magnons
In addition to the static magnetic properties described previously, the appearance of spin-wave excitations or so-called magnons is an important dynamic feature of magnetic materials, because their excitation spectrum can be compared directly with experimental results obtained, e.g., from neutron powder diffraction [161] or spin-resolved electron energy loss spectroscopy. [162,163] In collinear systems, the magnons are described by the excited states of the Heisenberg Hamiltonian (Equation (42)) under the assumption that the magnetic moments deviate only a little from their ground state directions. Their spectrum and the corresponding transverse magnetization deviations are given as the eigensystem [164] of the torque matrix which is associated with the Heisenberg Hamiltonian (Equation (42)) using also the magnetic spin moment S i and the exchange constants J ij . The e i are no vector components but e i ¼ þ1 when the moment in the ground state at site i points upwards or À1 otherwise. The index λ in Equation (46) labels eigenvalues. The e þ λi are interpreted as the spin-wave mode amplitude of mode λ at site i and ω λ is the spin-wave frequency space.
The calculation of the torque matrix in Equation (47) yields for many different materials a very good description of the magnon spectrum. For the TMO, e.g., the theoretical calculations compare well with the experimental results ( Figure 10). Some deviations can be attributed to the anisotropy and the alignment energy terms, which were neglected in our consideration. Furthermore, the theoretical curves generally underestimate the experimental energies in particular for larger wave vectors q away from the Γ point because of the underestimation of the exchange parameters, which directly contribute to the torque matrix (Equation (47)) (cf . Table 1). Moreover, this approach to calculate magnon spectra can be generalized also to disordered materials, which was demonstrated by Buczek et al. [163] for binary alloys.

Monte Carlo Simulations
In addition to the spin-wave spectra and the mean-field approximation of the magnetic transition temperature, the classical Heisenberg Hamiltonian (Equation (42)) can be used in a Monte Carlo (MC) simulation to estimate temperature-dependent quantities such as magnetic susceptibility, heat capacity, and magnetic phase transition temperature T m . Such a procedure is independent of the way of obtaining J ij .
We briefly summarize the important technical aspects but refer the reader to the book of Landau and Binder [165] for more details. With respect to the oxides, we simulated a lattice of the N magnetic sites repeating the systems unit cell L times in each direction. In general, the MC method can be applied to all kinds of systems from bulk, over thin films, to disordered lattices (see Section 6 for various examples on oxides). The latter might include random occupation of lattice sites with magnetic atoms. Periodic boundary conditions can be included and adapted for any geometry.
The desired observables are determined at each temperature step T from importance sampling as average over a large Figure 10. Calculated spin-wave dispersion of TMO obtained with HUTSEPOT (blue lines) and compared with experimental data (black diamonds and open circles). The paths in the BZ were chosen along several high-symmetry lines: from X ¼ ð0.25, 0.25, 0.25Þ via Γ ¼ ð0, 0, 0Þ to M ¼ ðÀ0.5, À 0.5, 0.5Þ, and then from M 0 via Γ to M 00 of the neighboring AFII BZ (see insets in the MnO panel). Adapted with permission. [103] Copyright 2009, American Physical Society.
www.advancedsciencenews.com www.pss-b.com number of Monte Carlo steps. Here, a single Monte Carlo step means in case of the Heisenberg Hamiltonian (Equation (42)) one sweep over all N sites in the lattice model. At each site i, a new orientation of its magnetic moment unit vector e 0 i is randomly determined. The variation of the total energy ΔE is calculated by evaluating Equation (42) for the old e i and the new orientation of the magnetic moment unit vector e 0 i . Afterwards, the acceptance of the new orientation is determined by means of the Metropolis algorithm. Then, the next site is considered and so on.
Before the observable is 'measured', we perform a large number of MC steps to reach the equilibrium state at a given temperature. After that, the average is taken over more MC steps. The optimal number of MC steps depends strongly on the type of considered system and its atomic species. Hence, the calculations are done several times with increasing number of MC steps to converge the results for the desired observables. The same has to be done for the lattice model size L due to finite-size effects or the number of included J ij from Equation (42).
We used here an implementation [64] of the MC method with the classical Heisenberg Hamiltonian specially adjusted for the use with HUTSEPOT. It was extended several times, [103,166,167] and applied to many different materials, especially oxides, as discussed in Section 6. It turns out that in particular for oxides and AFM materials additional observables can be defined, which can attribute for the correct T m . Here, we want to emphasize the staggered magnetization m s , as an average over the different AFM sublattices. [168] It was successfully applied to study the transition temperature in TMO [103] in combination with the fourorder cumulant U 4 . [169] Using the J ij obtained with the MFT for the TMO series (see Figure 9), the T m calculated with the MC method underestimate the experimental reference values ( Table 2). [103] The deviations being in the order of a few tens of kelvin fall in line with results for other oxide systems (see Section 6). Fischer et al. [103] argued with imperfect lattices in the experiments and not accounting for quantum nature or short-range order effects. Indeed, the theoretical determination of T m poses a complex problem trying to compare a perfect lattice of magnetic atoms with real macroscopic samples, which might have different grains or other possible lattice defects. In particular, we discuss subsequently how oxygen vacancies might increase T m in more complex oxides.

Selected Examples
We demonstrated above how the different aspects of our numerical methods work well for basic examples such as transition metal monoxides. The latter form rather simple lattice structures and many of their material properties are well understood. More complex oxide materials pose many more unresolved questions from the theoretical as well experimental point of view. Here, the proposed combination of theoretical methods proved to be very useful in the understanding of microscopic aspects sometimes not even accessible by experimental means: the magnetism in bulk ZnO [170,171] or at its surface, [172] the variation of magnetic properties with defects in the ferrite ZnFe 2 O 4 , [173,174] or the stability of hexagonal BaTiO 3 , [175] while ferroelectric BaTiO 3 can influence magnetic properties in on-top-deposited Co films. [115] In what follows, we focus on a few specific examples to demonstrate our multicode approach. Hereby, the most examples appeared either in perovskite or in double-perovskite structures with generic structure formula ABO 3 and ABB 0 O 6 , respectively. The perovskite prototype lattice model is cubic with B or B 0 cations at the center, an A cation at the corners, and oxygen at the faces of the cube, but strong internal distortions can appear depending mainly on the ratio of the cations A, B, and B 0 . We can compare, e.g., the lattice structure of the double perovskite Sr 2 FeMoO 6 (SFMO) with the heavily distorted perovskite structure of PrMnO 3 (PMO) (figures shown below). Such a lattice structure is also a general feature of other manganites, e.g., GdMnO 3 . [148,160] PMO was in particular interesting as the parent compound of the solid solution Pr 1Àx Ca x MnO 3 (PCMO). [177,178] 6.1. Magnetostructure of Pr 0.9 Ca 0.1 MnO 3 PCMO remains in the perovskite lattice structure over the whole concentration range 0 ≤ x ≤ 1 and reveals a series of different properties in dependence on x; i.e., it can be considered as a multifunctional semiconducting material. The magnetic structure for x ¼ 0.1 was characterized by Jirák et al. [179,180] by a canted AFM model. A change to x ¼ 0.2 (magnetic cluster glass behavior [181] ) or x ¼ 0.3 (insulator to metal transition [182,183] ) changes the situation completely. Thus, a critical reinvestigation of the magnetic structure for x ¼ 0.1 was performed by Tikkanen et al. [177] combining neutron powder diffraction, SQUID magnetometry, and first-principles calculations. The latter utilized two components out of our toolbox: 1) HUTSEPOT for the calculations of the electronic and magnetic properties and 2) the Monte Carlo method to determine temperature-dependent magnetic quantities. In addition, correlation corrections played an important role and were taken into account via the GGAþU method (Section 3.3). The mixing of Pr and Ca atoms in the cation sublattice was accomplished by means of the CPA (Section 4.1). More technical details of the calculations can be found in previous studies. [177,178] Figure 11 shows the DOS for x ¼ 0.1. It shows metallic behavior, which does not agree with experimental findings. Small concentrations of O vacancies are likely in the experimental samples. Assuming 2% oxygen vacancies restores the semiconducting character in the calculation.
Heisenberg exchange parameters can be effectively calculated based on the MFT, as discussed in Section 5.1.1. The results are shown in Figure 12. The coupling between Mn1 and Mn3, i.e., between two Mn planes, is given by J ¼ 4.8 meV, whereas the coupling between Mn1 and Mn4 is found to be AFM (J ¼ À1.0 meV). The effective coupling of two adjacent Mn planes is antiferromagnetic, as it is found by counting the interactions of the Mn atoms between the planes. MC simulations www.advancedsciencenews.com www.pss-b.com based on the Heisenberg exchange parameters lead to a Néel temperature of T N % 120 K. MC simulations with an applied magnetic field reveal that the Pr moments order simultaneously with the Mn moments, as found from experiments. The combination of neutron diffraction, SQUID magnetometry, and first-principles calculations leads to the conclusion that the lowtemperature magnetic moment of Pr is close to the high-spin value of 2μ B per formula unit, being 300% larger than previously thought.

Tuning of Magnetism in Sr 2 FeMoO 6 via O Vacancies
For the double perovskite SFMO, we were interested in similar properties as in the case of PCMO in the last section but with another motivation. SFMO shows in experimental and theoretical studies a high spin polarization. [184,185] This points to a half metallic character of SFMO, which provides a good starting point for spintronic devices, [7] and can be in fact well described by SIC calculations [186] or with a GGAþU approach. [187,188] We used the latter because DOS calculations with SIC by Szotek [186] show the d states of SFMO far below the Fermi energy, hinting at a too strong localization caused by the method. The spintronic devices would require SFMO in thin film form, which can be grown with high structural quality. [189] However, it is known that the preparation of SFMO thin films, e.g., with pulsed layer deposition, is sensitive to various kinds of defects (see for examples Figure 13). Hence, the magnetic transition temperatures T m of the thin films are found to be below those of bulk SFMO (T m % 400 K). [176] SFMO shows a ferrimagnetic ground state with Fe moments aligning antiparallel to the Mo moments.
In contrast to the Ca doping in PMO described earlier, SFMO shows already intrinsically some disorder at the two transition metal cation sites B and B 0 ; here Fe and Mo have almost equivalent lattice positions. Therefore, Mo might sit at the Fe site and vice versa-known as an antisite defect (ASD). It has been shown that those ASDs influence the electronic properties and reduce the spin polarization, [185,190] whereas T m was only affected for larger ASD concentrations. [191] We showed theoretically that more oxygen vacancies in SFMO will increase the magnetic transition temperature for the price of a reduced magnetic moment (Figure 14a). We combined again the magnetic exchange interactions J ij calculated with HUTSEPOT and the determination of T m by the Monte Carlo method. [167,187] The observed variation of magnetic moment and T m was attributed to the increasing AFM coupling between Fe and Mo, while at the same time the strength of the coupling between Fe─Fe pairs was reduced (Figure 14b). This observation is also supported by experiments with films additionally annealed after preparation. [176] The high temperatures in the annealing cause more oxygen vacancies, while in turn the measured T m is higher than for untreated samples.
In case of the thin films, the substrate will cause epitaxial strain in the SFMO films. In a joint experimental and theoretical study, we excluded epitaxial strain as a source varying T m on its own but found hints that the defect concentration varies substantially from the bulk sample to thin films. [176] Therefore, we studied with VASP the behavior of a number of point defects under a biaxial strain. [188] The internal structure of 2 Â 2 Â 1 supercells of (a) (b) Adapted with permission. [177] Copyright 2016, IOP Publishing Ltd.
(a) (c) (e) (f) Figure 13. Sketch of the SFMO structure as a 2 Â 2 Â 1 supercell: a,b) including internal relaxations calculated with VASP, c) oxygen vacancy, d) ASD, and e,f ) defect complexes. Adapted with permission. [176] Copyright 2016, American Chemical Society.   Figure 13). The total energies of the different supercells were then used to determine the defect formation energies as described by Nayak et al. [175] with where EðD, ε a Þ is the total energy of the supercell including an uncharged defect D. E host (ε a ) is the total energy of the host supercell of equal size under the influence of the strain ε a . Equation (48) has to be balanced by the chemical potentials μ i of the atomic type (i). The choice of μ i is related to the oxygenrich case (see the supporting information from Adeagbo et al. [188] ). In addition, the chemical potential for oxygen is connected to the oxygen partial pressure by thermodynamic relations. This allows a direct comparison of the theoretical results with the experimental growth conditions. The last term E corr in Equation (48) takes into account potential correction schemes for postprocessing formation energy results obtained by DFT methods. [175,193] Thus, we calculated the defect formation energy for different biaxial strains under oxygen-rich condition with Equation (48). The lattice constants in the xy of the supercell were varied and the resulting lattice constant in z direction relaxed. Those results were then related to a partial oxygen pressure corresponding to air (Figure 15a) and high vacuum (Figure 15b). At air partial pressure with a low E form ðV Sr , 0Þ, V Sr are the most likely defects, which is a little surprising, because the Sr stoichiometry is rarely taken into account as a potential source of variations in earlier studies. In contrast, E form ðV Sr , ε a Þ varies only a little for strains in an experimentally realizable range of AE3%. The formation energy for ASD increases for compressive strains, whereas it decreases slightly for tensile strains, resulting in an lower and higher probability of ASD, respectively ( Figure 15a). As E form ðASD, ε a Þ does not depend on the oxygen partial pressure, the curve remains at the same position also for lower oxygen partial pressure (Figure 15b).
Most importantly, oxygen vacancies are on the one hand very likely for low oxygen partial pressures. On the other hand, strain-independently from the direction-lowers the formation energy and, hence, increases the probability to form oxygen vacancies during the growing process. This will have in turn an impact on the magnetic properties with the relations discussed previously (Figure 14a).
In summary, we showed how sensitive the properties of SFMO are with different concentrations of defects. The latter can be directly influenced by epitaxial strain, demonstrating a path to tune the magnetic properties of oxides in general via variation of intrinsic point defects.

XMCD as a Chemical Fingerprint
Although we could show in the last section that oxygen vacanices in SFMO are one source of variations in the measured saturation magnetization, there could be many other sources which influence this single observable as well. Thus, we also study XAS and the XMCD. XAS and XMCD spectra are very sensitive to chemical composition, valency, and crystal structure and may, therefore, act as a fingerprint to the local chemical environment-most important in our studies on oxides. X-ray absorption is in general a many-body problem, but it can be approximated to some extent with in a one-particle picture. This can be better handled within first-principles methods, where it was formulated, e.g., in terms of the KKR Green's function method [58,194] or the LMTO method. [195,196] We apply the latter (see Section 2.2.2) to complement the study of electronic and magnetic properties within our multicode approach.
(a) (b) Figure 15. Defect formation energies obtained from VASP for different point defects in SFMO in dependence of biaxial strain at oxygen partial pressure equivalent to a) air and b) high vacuum. Adapted with permission. [188] Copyright 2018, American Physical Society.
(a) (b) Figure 14. Magnetic properties of oxygen-deficient Sr 2 FeMoO 6Àδ calculated with HUTSEPOT and CPA. a) Magnetic exchange interactions between nearest and next-nearest-neighbor Fe sites and between Mo and Fe sites. b) Total magnetic moment per f.u. and Curie temperature T C . Experimental saturation magnetization from Kircheisen and Töpfer. [192] Adapted with permission. [187] Copyright 2018, IOP Publishing Ltd.
www.advancedsciencenews.com www.pss-b.com In case of SFMO, the calculated XAS and XMCD spectra supported our investigation of oxygen vacancies with CPA. [115] A good agreement with experimental spectra at the L 2 and L 3 edge of Fe [197] was only possible when oxygen vacancies were also taken into account ( Figure 16). In fact, the small features at a and e in the experimental XAS results could not be explained without oxygen vacancies, which also cause the valency change of Fe from Fe 3þ to Fe 2þ . In adddition to SFMO, other double-perovskite materials were of interest, e.g., A 2 FeReO 6 with A ¼ Ca, Sr, and Ba. [198] Moreover, SFMO is not the only material discussed for spintronic applications. [7] Another promising candidate is the half-metallic La 0.7 Sr 0.3 MnO 3 , which is used as electrode in magnetic tunnel junctions. [6] Here, the interface to other oxide materials is crucial and we applied our multicode approach to study electronic and magnetic properties of SrRuO 3 /La 0.7 Sr 0.3 MnO 3 heterostructures. Although SrRuO 3 (SRO) can be correctly described from first-principles within the GGA approximation, [199] the Mn 3d electrons in La 0.7 Sr 0.3 MnO 3 (LSMO) are strongly localized and cannot be treated with the conventional DFT approach. At the same time, LSMO is a metal and the on-side Coulomb interaction can be efficiently screened by free electrons. To estimate the strength of correlation effects, we vary the Hubbard U eff (see Section 3.3) as a parameter to fit calculated XAS and XMCD spectra from the Mn L 2,3 edges with experimental observations (Figure 17). Due to restriction in the considered supercell, we approximate La 0.7 Sr 0.3 MnO 3 with La 0.75 Sr 0.25 MnO 3 . The best agreement with experiments is found for U eff ¼ 0.9 eV, whereas smaller (LDA) or larger (U eff ¼ 4.0 eV) values do not reproduce correctly the shape of the L 3 X-ray absorption spectrum (Figure 17a). Such a small U eff parameter also describes satisfactorily other properties of La 0.75 Sr 0.25 MnO 3 , e.g., the electronic structure or optical and magneto-optical Kerr (MOKE) spectra. [201] The XMCD spectrum at the Mn L 3 edge shows the prominent negative peak at 643.3 eV and a positive broad high-energy shoulder at 645-648 eV. Although all three theoretical results describe the major negative peak relatively well, U eff ¼ 0.9 eV covers the broad high-energy shoulder best. Thus, we conclude that the electronic correlations in the system are moderate and can be efficiently described with the relatively small Hubbard parameter U eff ¼ 0.9 eV, which was used in our works on the SRO/LSMO heterostructure. [202,203] Therein, SRO and LSMO are ferromagnetic but couple antiparallelly to each other over their interface (see Figure 18). The estimated strength of this AFM coupling is -13.3 meV and is caused by the spin polarization of the oxygen atoms at the interface, visible in the spin density. The latter is in SRO mostly of the same sign (blue color), whereas it alternates in sign along the Mn-O bonds in LSMO (see inset in Figure 18). In particular, at the interface, the oxygen ions are strongly polarized by the Mn atoms, leading to an asymmetry in the spin density. Induced by the negative spin density of the σ bond at the interface, the Ru atoms adjust their moments in parallel direction. Thus, we could show via the theoretical considerations why the SRO magnetization is oppositely aligned with respect to the LSMO magnetization.  Figure 16. Calculated a) left and b) right circularly polarized light XAS and c) XMCD results for the L 2,3 edge of Fe in SFMO with and without oxygen vacancies calculated with the LMTO method and compared to experimental measurements. [197] Adapted with permission. [187] Copyright 2018, IOP Publishing Ltd.  [200] www.advancedsciencenews.com www.pss-b.com

Microscopic Structure of CoO Bilayers on Ir(001)
The SFMO films discussed so far are "thin films" (a few 100 nm) from the experimental point of view, but the calculations were done for strained bulk material. In films of only a few monolayers, however, it is again different and the local structure can be even more important than before. The microscopic structure and electronic properties of CoO raised lengthy discussions across the literature. Bulk CoO is a wide-gap AFM insulator in the RS crystalline structure discussed earlier. CoO (111) is instead composed of alternating layers of Co and oxygen atoms (one bilayer) and should be actually polar. Hence, films of CoO (111) should be unstable despite their existence. This opened room for several speculations of charge compensation mechanisms and beyond. [204][205][206][207] It turned out that ultrathin CoO films compensate the additional charge by structural transition to a wurtzite-like structure, while the surface remains metallic. [204] In this study, we complemented surface X-ray diffraction (SXRD) and stress measurements by first-principles calculations with HUTSEPOT. The latter was especially suited for the task, as it is also designed for semi-infinite systems.
Two different films were grown by depositing Co atoms on an Ir(001) surface followed by annealing in an oxygen atmosphere. The structural characterization by SXRD showed different coverages of 1.6 and 2.0 bilayer of CoO (111). [204] All bond lengths are much smaller than for bulk CoO, whereas those inside the 2.0 ML film are a little bit smaller than those of the 1.6 ML film, to accommodate the larger coverage ( Figure 19). The other features of the two films are very similar and we only show here the important features for the 1.6 ML CoO film.
We adapted the experimental structure model to our first-principles simulations and applied the GGAþU method to describe the localized nature of Co d states to calculate electronic and magnetic properties. From the DOS we conclude that the films were metallic, which contributes to the compensation of the polarity ( Figure 20). Apparently, the thickness of the films is not large enough to open a bandgap. Our calculations indicate a very strong hybridization between cobalt d and oxygen p states throughout the whole valence band due to reduced distances between the Co and O atoms. This hybridization also enhanced the Co magnetic moments (see values in Figure 20).
The magnetic exchange interactions can be also easily calculated with Equation (43) for thin films. Thereby, not only the nearest and next-nearest J ij are nonnegligible as for CoO (Section 5.1), but many more are not small and have to be taken into account because the AFM super exchange coupling is significantly suppressed. The main exchange interactions have (a) (b) Figure 18. The MnO 2 /SrO-terminated interface in the La 0.7 Sr 0.3 MnO 3 / SrRuO 3 heterostructure. a) Structure model used as a periodic superlattice for calculations. Atomic types are color coded in the labels. b) Layerresolved magnetic moments obtained with HUTSEPOT. The inset shows the atomically resolved spin-density distribution with red/blue denoting spin up/down, respectively. Adapted with permission. [202] Copyright 2019, American Physical Society.
(a) (b) Figure 19. Lattice model for the 1.6 ML CoO film at Ir(001) derived from XRD measurements shown for a) top view and b) side view. The indices "i" and "s" mark atoms in the bilayer next to the Ir(001) substrate or at the surface, respectively. All bond lengths are given in angstroms. Adapted with permission. [204] Copyright 2014, American Physical Society.
www.advancedsciencenews.com www.pss-b.com AFM character, but many J ij are also FM or vanish completely due to the large variations in the Co─O bond lengths and the interaction with the underlying Ir substrate. Hence, we showed that the electronic and magnetic properties of the ultrathin CoO film are significantly altered, revealing a new character totally different from the wide-gap AFM CoO bulk. Those variations were strongly connected to the influence of the substrate and resulting structural variations within the ultrathin regime. Such deep understanding on the atomic level was hardly possible for experimental studies alone. In particular, the first-principles complemented the experimental findings nicely, which underlines their merits.

Conclusions
In this work, we presented a multicode approach for computational design of complex oxide materials which combines several complementary first-principles methods within the DFT. First, the crystalline structure can be obtained using pseudopotential codes, which are well designed for total energy calculations. Further information about the structure and chemical composition can be elucidated from first-principles simulations of observables, e.g., XAS and XMCD spectra, and their match to experimental results. The obtained structural information serve as input for calculations of the electronic and magnetic properties utilizing a first-principles Green's function method.
For the complex oxides, the strongly localized electrons are treated by means of the SIC-LSDA or LDAþU methods, while itinerant Bloch electronic states are described within the LDA or GGA approximations. Disorder effects can be treated efficiently using either a CPA based on a Green's function formalism or the SQS method, which mimics efficiently random alloy correlation functions up to next-nearest-neighbor interactions. The interplay of the different methods used in this multicode approach is illustrated by examples of electronic and magnetic structure studies of transition metal and complex oxides.
The approach is general and is successfully used to study electronic and magnetic properties of other classes of materials as well, e.g., diluted magnetic semiconductors, ferroelectrics, multiferroics, graphene-based compounds, topological insulators, Weyl semimetals, Heusler alloys, superconductors, metallic surfaces, interfaces, and clusters. Nevertheless, several applications are recently implemented to HUTSEPOT, in particular, the description of relativistic phenomena, or magnetic excitations in ordered and disordered materials. Further developments are currently being made to describe excited state properties beyond the conventional density functional applying the linear response theory [208,209] or the GW approximation. [210] [1] I.  The magnetic moments are shown for each Co atom. Adapted with permission. [204] Copyright 2014, American Physical Society.