Correlated materials design: prospects and challenges
CCorrelated Materials Design: Prospects and Challenges
Ran Adler , Chang-Jong Kang , Chuck-Hou Yee , and Gabriel Kotliar , Dept. of Physics & Astronomy, Rutgers, The State University of New Jersey, Piscataway, NJ 08854, USA Condensed Matter Physics and Materials Science Department,Brookhaven National Laboratory, Upton, New York 11973, USA (Dated: February 5, 2019)The design of correlated materials challenges researchers to combine the maturing, high through-put framework of DFT-based materials design with the rapidly-developing first-principles theory forcorrelated electron systems. We review the field of correlated materials, distinguishing two broadclasses of correlation effects, static and dynamics, and describe methodologies to take them intoaccount. We introduce a material design workflow, and illustrate it via examples in several materi-als classes, including superconductors, charge ordering materials and systems near an electronicallydriven metal to insulator transition, highlighting the interplay between theory and experiment witha view towards finding new materials. We review the statistical formulation of the errors of currentlyavailable methods to estimate formation energies. We formulate an approach for estimating a lower-bound for the probability of a new compound to form. Correlation effects have to be considered inall the material design steps. These include bridging between structure and property, obtaining thecorrect structure and predicting material stability. We introduce a post-processing strategy to takethem into account.
I. INTRODUCTION
The ability to design new materials with desired prop-erties is crucial to the development of new technology.The design of Silicon and Lithium-ion based materi-als are well known examples which led to the prolifer-ation of consumer hand-held devices today. Materialsdiscovery has historically proceeded via trial and error,with a mixture of serendipity and intuition being themost fruitful path. For example, all major classes ofsuperconductors–from elemental Mercury in 1911, to theheavy Fermions, Cuprates and most recently, the iron-based superconductors–have been discovered largely bychance .The dream of materials design is to create an effectiveworkflow for discovering new materials by combining ourtheories of electronic structure, chemistry and computa-tion. It is an inverse problem: start with the materialsproperties desired, and work to back out the chemicalcompositions and crystal structures which would lead todesirable properties. It requires a conceptual frameworkfor thinking about the physical properties of materials,and sufficiently accurate methods for computing them.In addition it requires algorithms for predicting crystalstructures and testing them for stability against decom-position, efficient codes implementing them and broadlyaccessible databases of known materials and their prop-erties.For weakly correlated materials, systems for whichband theory works, significant progress in all these frontshas been made. Fermi liquid theory justifies our think-ing of the excitations of a solid in terms of quasiparti-cles. Kohn Sham density functional theory (DFT) is agood tool for computing total energies and a good start-ing point for computing those quasiparticle propertiesin perturbation theory in the screened Coulomb interac-tions. Practical implementations of DFT such as LDA and GGA have become the underlying workhorse of thescientific community. Extensive benchmarks of softwareimplementations have shown that DFT reliably pro-duces the total energy of a given configuration of atoms,enabling comparisons of stability between different chem-ical polymorphs. The maturity of DFT, combined withsearchable repositories of experimental and calculateddata (Materials Project , OQMD , AFLOWlib andNIMS ), has fostered the growth of databases of com-puted materials properties to the point where one cansuccessfully design materials (see for example Refs. 7–9).Indeed these advances are beginning to pan out. Thesearch for new topological materials such as topologicalinsulators or Weyl semimetals is now greatly aided byelectronic structure calculations (for a recent review seeRef. 10). Another clear example of this coming of ageis the recent prediction of superconductivity in H S un-der high pressure near 190 K . Subsequently, hydrogensulfide was observed to superconduct near 200 K, thehighest temperature superconductor discovered so far .The situation is different for strongly correlated mate-rials. Many aspects of the physics of correlated electronmaterials are still not well understood. Correlated sys-tems exhibit novel phenomena not observed in weakly-correlated materials: metal-insulator transitions, mag-netic order and unconventional superconductivity aresalient examples. Designing and optimizing materialswith these properties would advance both technology andour understanding of the underlying physics.Furthermore, material specific predictive theory forthis class of materials is not fully developed, so even thedirect problem of predicting properties of correlated ma-terials with known atomic coordinates is very challeng-ing. It requires going beyond perturbative approaches,and we currently lack methods for reliably modeling ma-terials properties which scale up to the massive numberof calculations necessary for material design purposes. a r X i v : . [ c ond - m a t . s t r- e l ] F e b In this article, we examine the challenges of material de-sign projects involving strongly correlated electron sys-tems. Our goal is to present the state of the art in thefield stressing the outstanding challenges as it pertainsto correlated materials, and propose strategies to solvethem. We begin by providing a clear definition of cor-relations (Sec. II), distinguishing two important types,static and dynamic, and some available tools to treatthem. Next we introduce the material design workflow(Sec. IV). Then we give five examples of materials designin correlated systems to illustrate the application of ourideas (Sec. VI A-VI E) and conclude with a brief outlook.
II. WHAT ARE CORRELATED MATERIALS.STATIC AND DYNAMIC CORRELATIONS
The standard model of periodic solids views the elec-trons in a crystal as freely propagating waves with welldefined quantum numbers, crystal momentum and bandindex. Dating back to Sommerfeld and Bloch, it now hasa firm foundation based on the Fermi liquid theory andthe renormalization group, which explains why the ef-fects of Coulomb interactions disappear or “renormalizeaway” at low energies, and provides an exact descriptionof the excitation spectra in terms of quasiparticles. An-other route to the band theory of solid, is provided by thedensity functional theory in the Kohn-Sham implemen-tation, where a system of non-interacting quasiparticlesis designed so as to provide the exact density of a solid.While this wave picture of a solid has been extraordinar-ily successful and is the foundation for the descriptionof numerous materials, it fails dramatically for a class ofmaterials, which we will denote strongly correlated elec-tron systems.The basic feature of correlated materials is their elec-trons cannot be described as non-interacting particles.Since the constituent electrons are strongly coupled toone another, studying the behavior of individual parti-cles generally provides little insight into the macroscopicproperties of a correlated material. Often, correlated ma-terials arise when electrons are subjected to two compet-ing tendencies: the kinetic energy of hopping betweenatomic orbitals promotes band behavior, while the poten-tial energy of electron-electron repulsion prefers atomicbehavior. When a system is tuned so that the two energyscales are comparable, neither the itinerant nor atomicviewpoint is sufficient to capture the physics. The mostinteresting phases generally occur in this correlated anddifficult-to-describe regime, as we shall see in subsequentsections.These ideas have to be sharpened in order to quan-tify correlation strength, as there is no sharp boundarybetween weakly and strongly correlated materials. Ulti-mately one would like to have a methodology which canexplain the properties of any solid and which seeks tomake predictions for comparison with experimental ob-servations. To arrive at an operational definition of a correlated material, we examine DFT and how it relatesto the observed electronic spectra.The key idea behind DFT is that the free energy ofa solid can be expressed as a functional of the electrondensity ρ ( (cid:126)r ). Extremizing the free energy functional oneobtains the electronic density of the solid, and the valueof the functional at the extremum gives the total freeenergy of the material. The functional has the formΓ[ ρ ] = Γ univ [ ρ ] + (cid:82) d rV cryst ( r ) ρ ( r ) where Γ univ [ ρ ] is thesame for all materials, and the material-specific informa-tion is contained in the second term through the crys-talline potential. The universal functional is written as asum of T [ ρ ] the kinetic energy, E H , a Hartree Coulombenergy and a rest which is denoted as F xc the exchangecorrelation free energy. This term needs to be approx-imated since it is not exactly known, and the simplestapproximation is to use the free energy of the electrongas at a given density. This is called the Local DensityApproximation (LDA).The extremization of the functional was recast byKohn and Sham in the form of a single par-ticle Schr¨odinger equation with the Hartree atomicunits (cid:16) m e = e = (cid:126) = πε = 1 (cid:17)(cid:20) − ∇ + V KS ( (cid:126)r ) (cid:21) ψ (cid:126)kj ( (cid:126)r ) = (cid:15) (cid:126)kj ψ (cid:126)kj ( (cid:126)r ) . (1) (cid:88) (cid:126)kj | ψ (cid:126)kj ( (cid:126)r ) | f ( (cid:15) (cid:126)kj ) = ρ ( (cid:126)r ) (2)reproduces the density of the solid. It is useful to di-vide the Kohn-Sham potential into several parts: V KS = V H + V cryst + V xc , where one lumps into V xc exchangeand correlation effects beyond Hartree. In practice, theexchange-correlation term is difficult to capture, and isgenerally modeled by approximations known as the localdensity approximation (LDA) or generalized gradient ap-proximation (GGA). Density functional calculations us-ing the LDA/GGA approximation have become very pre-cise so that the uncertainties are almost entirely system-atic. To get a feel for the numbers, convergence criteriaof 10 − to 10 − meV/atom are routinely used whereasdifferences between experimental and theoretical heatsof formation routinely differ by over 100 meV/atom .The eigenvalues (cid:15) (cid:126)kj of the solution of the self-consistentset of Eqs. (1) and (2) are not to be interpreted as exci-tation energies. Konh Sham excitations are not Landauquasiparticle excitations. The latter represent the exci-tation spectra, which are the experimental observable inangle resolved photoemission and inverse photoemissionexperiments and should be extracted from the poles ofthe one particle Green’s function: G (cid:16) (cid:126)r, (cid:126)r (cid:48) , ω (cid:17) = 1 ω + ∇ + µ − V H − V cryst − Σ (cid:16) (cid:126)r, (cid:126)r (cid:48) ,ω (cid:17) . (3) FIG. 1. Schematic diagrams for the GW method. Startingfrom some G a polarization bubble is constructed, which isused to screen the Coulomb interactions resulting in an inter-action W. This W is then used to compute a self-energy Σ GW using W and G . To obtain the full Green’s function G inEq. (3), one goes from Σ GW to Σ by subtracting the nec-essary single particle potential and uses the Dyson equation G − = G − − Σ as discussed in the text. Adapted from 15.
Here µ is the chemical potential and we have singled outin Eq. (3) the Hartree potential V H ( (cid:126)r ) = (cid:82) ρ ( (cid:126)r (cid:48) ) | (cid:126)r − (cid:126)r (cid:48) | d r (cid:48) ex-pressed in terms of the exact density, V cryst is the crystalpotential and we lumped the rest of the effects of thecorrelation in the self-energy operator Σ (cid:16) (cid:126)r, (cid:126)r (cid:48) ,ω (cid:17) whichdepends on frequency as well as on two space variables.Since taking Σ = V xc generates the Kohn-Sham spec-tra, we define weakly correlated materials as ones where | Σ( ω ) − V xc | (4)is small for low energies, so our definition of weakly corre-lated materials are those for which the Kohn-Sham spec-tra is sufficiently close to experimental results.We can refine this definition, by taking into accountfirst order perturbation theory in the screened Coulombinteractions, taking LDA as a starting point. This is the G W method, which we now describe using diagrams inFig. 1. This figure first describes the evaluation of thepolarization bubble ΠΠ ( t, t (cid:48) ) = G ( t, t (cid:48) ) G ( t (cid:48) , t ) (5)Next, the screened Coulomb potential W in the randomphase approximation (RPA) which is the infinite sum ofdiagrams depicted and represent the expression W − = v − Coul − Π (6)where v Coul is the bare Coulomb potential. Then oneproceeds to the evaluation of a self-energyΣ GW = G W (7) which represents the lowest order contribution in pertur-bation theory in W (given in real space by Fig. 1), andthen G − = G − − Σ using Dyson’s equation. G above is just a Green’s function of non-interactingparticles, and it can thus be defined in various ways, lead-ing to different variants of the GW method. In the “one-shot” (that is, a method with no self-consistency loop)GW method (aka G W ) one uses the LDA Kohn-ShamGreen’s function. G ( iω ) − = iω + µ + 12 ∇ − V H − V cryst − V LDAxc . (8)and the self-energy is thus taken to be Σ = Σ GW − V LDAxc .Through this paper we use a matrix notation loosely andview operators as matrices. For example, in the Dysonequations
W, v
Coul , Π are operators (matrices) with ma-trix elements (cid:104) r | W ( ω ) | r (cid:48) (cid:105) = W ( ω, r, r (cid:48) ), (cid:104) r | Π( ω ) | r (cid:48) (cid:105) =Π( ω, r, r (cid:48) ), (cid:104) r | v Coul | r (cid:48) (cid:105) = | r − r (cid:48) | .This G W method , introduced by Hybertsen andLouie , systematically improves the gaps of all semicon-ducting materials. We show this in Fig. 2. The successof this G W method implies that in this kind of materi-als the Kohn-Sham references system is sufficiently closeto the exact self-energy that the first order perturbationtheory correction Σ G W ( ω ) − V LDA/GGA xc brings us closeenough to the experimental results.However, there are many materials (usually containingatoms with open d or f shells) where the photoemissionspectra (and many other physical properties) are not sowell described by this method. A successful many bodytheory of the solid state aims to describe all these sys-tems. For the most widely used DFT starting points,LDA and GGA, what is the physical basis for the devia-tions in Eq. (4)? It is useful to think about two limitingcases, one in which the self-energy Σ is a strong func-tion of frequency, in which case we talk about dynami-cal effects, and a case where Σ is strongly momentum-dependent, or in real space - highly non-local, but weaklyfrequency dependent, and we talk about static correla-tions. In materials with strong dynamical correlationsthe spectral function displays additional peaks, which arenot present in the band theory, and reflect the atomicmultiplets of the material.Electron correlation is customarily divided into dy-namical and non-dynamical, but there is no strict defini-tion of these terms. In the context of Quantum Chem-istry calculations, these terms are mainly used to describethe ability of different methods to capture significantcorrelation effects, and the type of wave function whichwould approximate the exact solution of the Schr¨odingerequation. Non-dynamical or static correlations in thechemistry context means that energetically-close / degen-erate electronic configurations are appreciably present inthe wave function. This requires multiple Slater deter-minants of low lying configurations, and multi-referencemethods to describe them, such as the multi-referenceHartree Fock method, or multi-reference coupled clus-ter methods. Dynamical correlation refers to a situation FIG. 2. Theoretically-determined semiconductor gap in a one shot LDA G W calculation versus experiment (data compliedby E. Shirley). Adapted from Chapter III. “First-principles theory of electron excitation energies in solids, surface, and defects”(article author: Steven G. Louie) in Topics in Computational Materials Science, edited by C. Y. Fong (World Scientific, 1998)[17]. Diamonds are the G W excitation gap, while the crosses are the LDA value. where a single Slater determinant, such as a closed shellconfiguration of some orbitals, is a good reference system- which then needs to be dressed by including double (orhigher) excitations from strongly occupied core shells toempty orbitals. In addition, other virtual processes canmodify the orbitals of the original slater determinant.This situation is well described in the standard coupledcluster method which is considered the gold standard inQuantum Chemistry .Confusingly, the chemist’s delocalization error corre-sponds to our definition of k dependent self-energy, whichwe denote as static correlation (since it does not in-volve frequency dependence of the self-energy), while thechemist’s static correlation corresponds to what we calldynamical correlation as it requires a strongly frequencydependent self-energy in condensed matter physics. Weuse the solid state physicist convention in this article.Another useful way to classify the correlations is by thelevel of locality of the self-energy. Introducing a complete basis set of localized wave functions labeled by site andorbital index we can expand the self-energy asΣ ( (cid:126)r, (cid:126)r (cid:48) , ω ) = (cid:88) α (cid:126)R,β (cid:126)R (cid:48) χ ∗ α (cid:126)R ( (cid:126)r )Σ α (cid:126)R,β (cid:126)R (cid:48) ( ω ) χ β (cid:126)R (cid:48) ( (cid:126)r (cid:48) ) . (9)The self-energy is approximately local when the on-site term R = R (cid:48) in Eq. (9) is much larger than the rest.Notice that the notion of locality is defined with referenceto a basis set of orbitals.Equation (9) allows us to introduce an approximationto the self-energy involving a sum of a non-local butfrequency independent term plus a frequency dependentbut local self-energy:Σ( (cid:126)k, ω ) (cid:39) Σ( (cid:126)k ) + (cid:88) (cid:126)R,αβ ∈ L | (cid:126)Rα (cid:105) Σ α (cid:126)R,β (cid:126)R ( ω ) (cid:104) (cid:126)Rβ | (10)This ansatz was first introduced by Sadovskii et al . It isuseful when the sum over orbitals in Eq. (10) runs over a O r de r O f I n t e r a c t i on (non) locality B a s i s S e t s i z e S i ng l e - s i t e D M FT - c l u s t e r D M FT FIG. 3. Two complementary approaches to the treatment ofcorrelations. One axis represents the systematic perturba-tive expansion in powers of an interaction (for example thescreened Coulomb interaction W). The second axis, sums per-turbation theory to all orders, but at the local level. Whenlocality is just a lattice site, we have the single site DMFT,improvements involve larger clusters. In addition when onegoes beyond model Hamiltonians towards the realistic treat-ment of solids we need to introduce a basis set and estimatethat the results are converged as function of the size of thebasis set. small set L (much smaller than the size of the basis set),for example over a single shell of d or f orbitals. Thisform captures both static and dynamical correlations andis also amenable to computation using Dynamical MeanField Methods to be introduced in section III. III. HOW TO TREAT CORRELATIONS
Having defined correlations as a departure of theGreen’s function from the results of lowest order per-turbation theory around LDA (i.e. G W ), we now re-view various ways to correlations into account. Oneshould keep in mind that different materials may requirestronger momentum or frequency dependence in theirself-energy, and may exhibit different degrees of locality.This section lays out several complementary approachesto treat correlations beyond G W . They represent dif-ferent compromises between speed and accuracy, and cantarget different levels of locality and different correlationsstrengths. A schematic view of the grand challenge posedby the treatment of correlations in the solid state is pre-sented in Fig. 3, which explains the need to converge thecalculations along multiple axis. Linearized Self Consistent Quasiparticle GW .We begin our treatment with the GW approximation,which was introduced in the previous section. One obvi-ous flaw of the G W method is its dependence on theLDA input. This makes the method increasingly inaccu-rate as the strength of the correlations increase. One wayto eliminate this dependence, is to introduce some level of self consistency. Hedin proposed a full self-consistent GW scheme, namely to use G = G is in Eq. (7). We canthink of this as setting V xc = 0, so it is not used in inter-mediate steps. There are numerous advantages, however,to using a non-interacting form for G in the algorithm,and in practice the spectra in self-consistent GW turnedout to be consistently worse in solids than the non-self-consistent approach for spectral properties . Neverthe-less, GW can be reasonably accurate for total energy cal-culations, as they can be obtained as stationary pointsof a functional .To improve on the spectra relative to G W while re-taining some level of self consistency so as not to dependon the starting point, the self-consistent quasi-particleGW (QPGW) was proposed. Here one uses the “best”non-interacting Green’s function G , which is defined interms of an “exchange and correlation potential” V QP GWxc chosen to reproduce the spectra of the full G as closelyas possible: G QP GW ( iω ) − = iω + µ + 12 ∇ − V H − V cryst − V QP GWxc . (11)To determine V QP GWxc (which once again we viewas a matrix with matrix elements (cid:104) r | V QP GWxc | r (cid:48) (cid:105) = V QP GWxc ( r ) δ ( r − r (cid:48) )), it was proposed to approximate thespectra and the eigenvectors of G by those of G QP GW - bysolving a set of non-linear equations on the real axis .An alternative approach that works on the imaginaryaxis is to linearize the GW self-energy at each iteration.Namely, after the evaluation of the self-energy in Eq. (7),this quantity is Taylor expanded around zero frequency(hence the name “linearized”):Σ lin ( (cid:126)k, iω ) = iω (1 − Z ( (cid:126)k ) − ) + Σ( (cid:126)k, G QP GW ( iω ) is obtained by solving the usual Dysonequation with the linearized self-energy, and multiplyingthe result by the quasiparticle residue, Z , to obtain aproperly normalized quasiparticle Green’s function. G QPGW0 = (12) (cid:113) Z ( (cid:126)k )[ iω + µ + 12 ∇ − V H − V cryst − Σ lin ] − (cid:113) Z ( (cid:126)k )Note that this defines the exchange correlation poten-tial of the self-consistent QPGW method. This method,the linearized self-consistent quasiparticle GW, was intro-duced in Ref. 26 and an open source code to implementthis type of calculation in the linearized augmented planewave (LAPW) basis set is available in Ref. 27.The GW or RPA method captures an important phys-ical effect. Electrons are charged objects which interactvia the long range Coulomb interactions. Quasiparti-cles, on the other hand, interact through the screenedCoulomb interaction. They are composed of electronssurrounded by screening charges, thus reducing the FIG. 4. Dynamical Mean Field Theory (DMFT) maps(or truncates ) a lattice model to a single site embeddedin a medium (impurity model) with a hybridization strengthwhich is determined self consistently. Adapted from Ref. 28. strength and the range of their interaction. For this rea-son, in many model Hamiltonians describing metals, onlythe short range repulsion is kept. On the other hand, it iswell known that the RPA fails in describing the pair cor-relation function at short distances. One can say that theGW method captures the long range of the screening ef-fects of the long range Coulomb interactions and producea self-energy which is non-local in space, but with a weakfrequency dependence (indeed the self-energy is linear ina broad range of energies). It turns out that this methodis not able to capture the effects of the short range partof the Coulomb interactions which in turn induces strongfrequency dependence (i.e. strong non-locality in time),but in turn is much more local in space.
Dynamical Mean Field Theory (DMFT) . To cap-ture dynamical local correlations one uses Dynamicalmean field theory , which is the natural extension of theWeiss mean field theory of spin systems to treat quantummechanical model Hamiltonians. Dynamical Mean FieldTheory becomes exact in the limit of infinite dimensions,which was introduced by Metzner and Vollhardt . Withsuitable extensions it plays an important role in realis-tically describing strongly correlated electron materials.Here we describe the main intuitive DMFT ideas as aquantum embedding, starting from the example of a one-band Hubbard model (describing s electrons), in whichthe relevant atomic configurations are | (cid:105) , |↑(cid:105) , |↓(cid:105) , |↑↓(cid:105) asdescribed in Fig. 4. It involves two steps. The first step,focuses on a single lattice site, and describes the rest ofthe sites by a medium with which an electron at this sitehybridizes. This truncation to a single site problem iscommon to all mean field theories. In the Weiss meanfield theory one selects a spin at a given site, and re-places the rest of the lattice by an effective magnetic fieldor Weiss field. In the dynamical mean field theory, thelocal Hamiltonian at a given site is kept, and the kineticenergy is replaced by a hybridization term with a bath ofnon-interacting electrons, which allows the atom at theselected site to change its configuration. This is depictedin Fig. 4 where we apply the method to the one-bandHubbard model. The system consist of one band of s FIG. 5. The DMFT impurity model is used to generate irre-ducible quantities such as self-energies and one particle ver-tices. These are then embedded in the lattice model to gen-erate momentum dependent lattice quantities such as spectralfunctions, or spin susceptibilities. Adapted from 15. electrons. The Fourier transform of the hopping integralis given by t ( −→ k ).It is used in the second step, which involves the re-construction of lattice observables by embedding thelocal impurity self-energy into a correlation function ofthe lattice, G latt ( (cid:126)k, iω ) − = iω + µ − t ( (cid:126)k ) − Σ imp ( iω ) . Here Σ imp ( iω ) are viewed as functionals of the Weissfield. The requirement that (cid:80) k G latt = G loc determinesthe Weiss field. Table I summarizes the analogies be-tween Weiss mean field theory and dynamical mean fieldtheory. Weiss Mean Field Theory Dynamical Mean Field TheoryIsing Model → Single Spin Hubbard Model → in effective Weiss Field Impurity in effective bathWeiss field: h eff effective bath: ∆( ıω n )Local observable: m = < s i > Local Observable: G loc ( ıω n )Self-consistent condition: Self-consistent condition:tanh (cid:16) β (cid:80) j J ij s j (cid:17) = m iω n − E imp − ∆ ( iω n ) − Σ ( iω n ) = (cid:2)(cid:80) (cid:126)k G (cid:126)k ( iω n ) (cid:3) − TABLE I. Corresponding quantities in Dynamical MFT(right) and Weiss or static MFT in statistical mechanics (left).
The DMFT mapping of a lattice model into an impu-rity model gives a local picture of the solid, in terms ofan impurity model, which can then be used to generatelattice quantities such as the electron Green’s functionand the magnetic susceptibility by computing the cor-responding irreducible quantities. This is illustrated inFig. 5.The self-consistent loop of DMFT is summarized in thefollowing iterative cycle E imp , ∆ ( iω n ) → Impurity Solver → Σ imp ( iω n ) , G loc ( iω n ) ↑ ↓ G (cid:126)k ( iω n ) =Truncation ← iωn + µ − t ( (cid:126)k ) − Σ( iωn ) ← Embedding
From the simplest model, the one-band Hubbardmodel, one can proceed to more realistic descriptions ofcorrelated materials by replacing t ( (cid:126)k ) by a tight-bindingmodel Hamiltonian matrix. The DMFT equations canbe derived from a functionalΓ DMF T model (cid:104) G αβ, (cid:126)R , Σ αβ, (cid:126)R (cid:105) = (13) − Tr ln (cid:104) iω n − H ( (cid:126)k ) − Σ αβ (cid:126)R ( iω n ) (cid:105) − (cid:88) n T r [Σ ( iω n ) G ( iω n )] + (cid:88) (cid:126)R Φ (cid:104) G αβ, (cid:126)R , U (cid:105) where Φ (cid:104) G αβ, (cid:126)R , U (cid:105) is the Baym-Kadanoff functional -the sum of all two particle irreducible diagrams in termsof the full Green’s function G , and the Hubbard inter-action U , which denotes a four rank tensor U αβγδ . Itcan also be evaluated from the Anderson impurity modelexpressed in terms of the full local Green’s function ofthe impurity G . The impurity model is the engine of aDMFT calculation. Multiple approaches have been usedfor its solution, and full reviews have been written onthe topic. The introduction of continuous time MonteCarlo method for impurity models have provided nu-merically exact solutions reducing the computational costrelative to the Hirsch-Fye algorithm that was used in ear-lier DMFT studies. DFT+DMFT method.
This is the next step to-wards a more realistic description of solids. It was in-troduced in Refs. 32 and 33. In these early implemen-tations, it consisted of replacing the Hamiltonian H ( (cid:126)k )by the Kohn-Sham matrix in Eq. (13) with a correctionto subtract the correlation energy that is contained inthe Kohn-Sham Hamiltonian (double counting correc-tion). The original DFT calculations were carried outwith an LDA exchange and correlation potential, butthey could be done with GGA and other functionals. Fur-thermore the exchange and correlation potential in theDyson equation for the Green’s function can be replacedby another static mean field theory like hybrid DFT orQPGW, but in the following we will use the terminologyLDA+DMFT.Starting from the Anderson model Hamiltonian pointof view, one divides the orbitals into two sets. The firstset contains the large majority of the electrons are prop- erly described by the LDA Kohn-Sham matrix. The sec-ond set contains the more localized orbitals ( d -electronsin transition metals and f -electrons in rare earths andactinides) which require the addition of DMFT correc-tions. A subtraction (called the double counting correc-tions) takes into account that the Hartree and exchangecorrelation has been included in that orbital twice, sinceit was treated both in LDA and in DMFT. The earlyLDA+DMFT calculations, proceeded in two steps (one-shot LDA+DMFT). First an LDA calculation was per-formed for a given material. Then a model Hamiltonianwas constructed from the resulting Kohn-Sham matrixcorrected by E dc written in a localized basis set. The val-ues of a Coulomb matrix for the correlated orbitals wereestimated or used as fitting parameters. Finally DMFTcalculation were performed to improve on the one particleGreen’s function of the solid.In reality, the charge is also corrected by the DMFTself-energy, which in turn changes the exchange and cor-relation potential away from its LDA value. Thereforecharge self-consistent LDA+DMFT is needed. This wasfirst implemented in Refs. 34 and 35.For this purpose it is useful to notice that theLDA+DMFT equations can be derived as stationarypoints of an LDA+DMFT functional, which can beviewed as a functional of the density and local Green’sfunction of correlated orbitals. This is a spectral densityfunctional. Evaluating the functional at the stationarypoint gives the free energy of the solid, and the station-ary Green’s functions gives us the spectral function of thematerial. We can arrive at the DFT+DMFT functionalby performing the substitutions − ∇ + V KS ( (cid:126)r ) for H ( (cid:126)k )in the model DMFT functional Eq. (13) and then addingterms arising from the density functional theory, namely:Γ DF T + DMF T (cid:104) ρ ( (cid:126)r ) , G αβ, (cid:126)R , V KS ( (cid:126)r ) , Σ αβ, (cid:126)R (cid:105) = Γ DMF T model [ H ( (cid:126)k ) → − ∇ + V KS ( (cid:126)r )]+ Γ [ V KS ( (cid:126)r ) , ρ ( (cid:126)r )] − Φ DC (14)whereΓ [ V KS ( (cid:126)r ) , ρ ( (cid:126)r )] = − (cid:90) V KS ( (cid:126)r ) ρ ( (cid:126)r ) d r + (cid:90) V ext ( (cid:126)r ) ρ ( (cid:126)r ) d r + 12 (cid:90) ρ ( (cid:126)r ) ρ ( (cid:126)r (cid:48) ) | (cid:126)r − (cid:126)r (cid:48) | d rd r (cid:48) + E DF Txc [ ρ ](15)We then arrive at the DFT+DMFT functional whichwe write in full below.Γ DF T + DMF T (cid:104) ρ ( (cid:126)r ) , G αβ, (cid:126)R , V KS ( (cid:126)r ) , Σ αβ, (cid:126)R (cid:105) = − Tr ln iω n + µ + ∇ − V KS − (cid:88) R,αβ ∈ L χ ∗ α (cid:126)R ( (cid:126)r ) Σ αβ (cid:126)R ( iω n ) χ β (cid:126)R ( (cid:126)r (cid:48) ) − (cid:90) V KS ( (cid:126)r ) ρ ( (cid:126)r ) d r − (cid:88) n Tr [Σ ( iω n ) G ( iω n )] + (cid:90) V cryst ( (cid:126)r ) ρ ( (cid:126)r ) d r + 12 (cid:90) ρ ( (cid:126)r ) ρ ( (cid:126)r (cid:48) ) | (cid:126)r − (cid:126)r (cid:48) | d rd r (cid:48) + E DF Txc [ ρ ] + (cid:88) (cid:126)R Φ (cid:104) G αβ, (cid:126)R , U (cid:105) − Φ DC . (16)Φ is the sum of two-particle irreducible diagrams writ-ten in terms of G and U . It was written down first inRef. 34, building on the earlier work of Chitra and Kotliar . It is essential for total energy calculations whichrequire the implementation of charge self-consistency inthe LDA+DMFT method. The first implementation ofcharge self-consistent LDA+DMFT was carried out in afull potential linear muffin-tin orbital (FP-LMTO) basisset . It was used to compute total energy and phononsof δ -plutonium .Alternatively, one can include the hybridization func-tion ∆ or the Weiss field G as an independent variable inthe functional in order to see explicitly the free energyof the Anderson Impurity Model, G − αβ, −→ R = G − atom α,β, −→ R − ∆ αβ, −→ R : F imp (cid:104) G − αβ, (cid:126)R (cid:105) = − ln (cid:90) D [ c † c ] e − S imp [ c † ,c ] with S imp [ G − αβ, (cid:126)R ] = − (cid:88) αβ (cid:90) dτ dτ (cid:48) c † α ( τ ) G − αβ, (cid:126)R ( τ, τ (cid:48) ) c β ( τ (cid:48) )+ U αβγδ (cid:90) dτ c † α ( τ ) c † β ( τ ) c δ ( τ ) c γ ( τ )So that:Γ DF T + DMF T (cid:104) ρ ( (cid:126)r ) , G αβ, (cid:126)R , V KS ( (cid:126)r ) , Σ αβ, (cid:126)R , αβ, (cid:126)R , G αβ, (cid:126)R (cid:105) = − Tr ln iω n + µ + ∇ − V KS − (cid:88) R,αβ ∈ L χ ∗ α (cid:126)R ( (cid:126)r ) Σ αβ (cid:126)R ( iω n ) χ β (cid:126)R ( (cid:126)r (cid:48) ) − (cid:90) V KS ( (cid:126)r ) ρ ( (cid:126)r ) d r − (cid:88) Tr (cid:2) ( G − − Σ ( iω n ) − G − ) G ( iω n ) (cid:3) + (cid:90) V cryst ( (cid:126)r ) ρ ( (cid:126)r ) d r + 12 (cid:90) ρ ( (cid:126)r ) ρ ( (cid:126)r (cid:48) ) | (cid:126)r − (cid:126)r (cid:48) | d rd r (cid:48) + E DF Txc [ ρ ] + (cid:88) (cid:126)R F imp (cid:104) G − αβ, (cid:126)R (cid:105) − Tr ln[ G αβ, (cid:126)R ] − Φ DC . (17)The form of the LDA+DMFT functional makes itclear that the method is independent of the basis setused to implement the electronic structure calculation,provided that the basis is complete enough. On theother hand, it is clearly dependent on the parameter U chosen, on the form of the double counting correc-tion and the choice of the projector (i.e., the orbitals χ α ( (cid:126)r ) with α ∈ L that enter this definition) and theexchange correlation functional E DF Txc . A projector ofthe form P ( r, r (cid:48) ) = (cid:80) αβ ∈ L χ ∗ α −→ R ( −→ r ) χ β −→ R ( −→ r (cid:48) ) was usedto define a truncation from G to G loc . The inverse of P is the embedding operator E defined by P · E = I L where I L is the identity operator in the space spannedby the correlated orbitals. If one restricts E · P to thespace L , one also obtains the identity operator in thatspace. E is used to define an embedding of the self-energyΣ( r, r (cid:48) ) = (cid:80) αβ E α,β ( r, r (cid:48) )Σ locα,β . However, more general projectors can be consideredas long as causality of the DMFT equations is satisfied.Ideas for choosing an optimal projector for LDA+DMFTbased on orbitals were presented in Ref. 39. Choosingsuitable projectors (and correspondingly a suitable valueof the U matrix and a proper double counting correc-tion) is crucial for the accuracy of an LDA+DMFT cal-culation as demonstrated recently in the context of thehydrogen molecule . DFT+DMFT is now a widely usedmethod. It has been successfully used across the periodictable, and has been implemented in numerous codes .Still there is ample room for advances in implementation,and on providing a firm foundation of the method. Onecan view the DFT+DMFT functional written above, asan approximation to an exact DFT+DMFT functional,which would yield the exact density and spectra of thesolid . This viewpoint has been used recently to provide FIG. 6. Hedin’s Equations give an exact representation of thecorrelation function of the Bosonic and Fermionic correlationin an expansion in G and W. They can be obtained by settingthe functional derivatives of the Eq. (18) to 0. Φ(
G, W )(first line) is the set of 2-PI skeleton diagrams (in G andW), where by convention the symmetry weights are omitted.The derivative by G (line 2) shows how the self-energy Σ isdefined in terms of a 3-legged vertex Γ. The derivative by W(line 3) equals the polarization Π . The bottom line shows thedefinition of the vertex Γ . an expression for the double counting correction Φ DC .An alternative perspective goes back to a fully diagram-matic many body theory of the solid, and examines howDFT+DMFT would fit in that framework as an approx-imation. We turn to this formulation next. Fully Diagrammatic Methods:
The free energyof the solid can also be expressed as a functional of G ( x, x (cid:48) ) and W ( x, x (cid:48) ) by means of a Legendre trans-formation and results in Refs. 36 and 52, where E H = (cid:82) ρ ( −→ r ) ρ ( −→ r (cid:48) ) | −→ r −−→ r (cid:48) | d rd r (cid:48) , Φ is defined as sum of all 2-particleirreducible diagrams which cannot be divided into twoparts by cutting two Green’s functions lines (which canbe either G’s or W’s):Γ [ G, W, Σ , Π] = − Tr ln (cid:2) G − − Σ (cid:3) − Tr [Σ G ]+ 12 Tr ln (cid:2) v − Coul − Π (cid:3) −
12 Tr [Π W ]+ E H + Φ [ G, W ] , (18)This reformulation is exact and leads to the exact Hedin’sequations, shown in Fig. 6. To convert this generalmethod into a tool of practical use, strong approxima-tions have to be introduced.The lowest order graphs of Eq. (18), shown in Fig. 7,reproduce the self-consistent GW approximation: takingfunctional derivatives of the low order functional withrespect to the arguments produces the same equations asthe GW approximation.To summarize the discussion so far, we recall that forsemiconductors, non-local (but weakly frequency depen-dent) correlation effects are needed to increase the gapfrom its LDA value. This admixture of exchange, canbe done within the GW method, or using hybrid densityfunctionals. It reflects the importance of exchange be-yond the LDA, which is due to the long-range but static FIG. 7. Lowest order graphs in the Φ-functional of Eq. (18).They give rise to the fully self-consistent GW approximation,as the saddle point equations. Note that the first term hereis the Hartree energy. Adapted from 15. part of the Coulomb interaction. These are the staticcorrelation effects. It has recently been shown, thatthis type of correlation effects are important in materi-als near a metal-to-insulator transition such as BaBiO or HfClN and can have a dramatic effect in enhancingthe electron-phonon interaction relative to its LDA esti-mated value. In these systems, a strongly k-dependentself-energy effect, Σ( k ), is much more important thanfrequency dependence, and here GW methods work well.On the other hand frequency dependence, and its im-plied non-locality in time, is crucial to capture Mott orHund’s physics. This physics tends to be local in spaceand can be captured by DMFT. Static mean field theo-ries, such as the LDA, do not capture this non-locality intime, and therefore fail to describe Mott or Hund’s phe-nomena. DFT+DMFT can treat strong frequency de-pendency, but has k-dependence only as inherited fromthe k-dependence of DFT exchange and correlation po-tential, the k-dependence of the embedding and the dou-ble counting shift.In real materials both effects are present to some de-gree - thus motivating physically the ansatz, Eq. (10).Some examples discussed recently are Ce O (using hy-brid DFT+DMFT) in Ref. 54 and the iron pnictides andchalcogenides in Ref. 19.We now describe a route proposed by Chitra toembed DMFT into a many-body approach to electronicstructure within a purely diagrammatic approach formu-lated in the continuum.If one selects a projector, which allows us to definea local Green’s function, it was suggested in Refs. 36,37, and 55 that one can perform a local approximationand keep only the local higher order graphs in a selectedorbital.Φ [ G, W ] (cid:39) Φ EDMF T [ G loc , W loc , G nonlocal = 0 , W nonlocal = 0] +Φ GW − Φ GW [ G loc , W loc , G nonlocal = 0 , W nonlocal = 0]Since the lowest-order graph is contained in the GWapproximation, one should start from the second ordergraph and higher order. This Φ GW + DMF T functional isshown in Fig. 8.These ideas were formulated and fully implemented inthe context of a simple extended Hubbard model .0 FIG. 8. Comparison of the functionals for the methods de-scribed in the text. The Hartree diagram was dropped sinceit appears in all methods . An open problem in this area, explored in Ref. 58 is thelevel of self consistency that should be imposed. This im-portant issue is already present in the implementation ofthe GW method, and the work of Ref. 58 should be revis-ited using the lessons from the QPGW method . Therehas been a large number of works exploring GW+DMFTand related extensions and combinations, and we referthe reader to recent reviews for the most recent refer-ences . Recently, we proposed the self consistent quasi-particle GW+DMFT , as a theory that contains boththe most successful form of a GW approximation andDMFT as limiting cases. Further understanding of thismethod requires the systematic treatment of vertex cor-rections, an approach which is now vigorously pursued.The LDA+U is a method was introduced by Anisi-mov et al. . It was made rotationally invariant inRefs. 63 and 64. One can view this is as a special case ofLDA+DMFT, where in the functional (Eq. (16)) Φ (thesum of graphs) is restricted to the Hartree-Fock graphs:Φ → Φ HF , and Σ( ıω n ) is replaced by a constant matrix λ . Then the LDA+U functional Γ LDA + U can be writtenin as follows:Γ LDA + U [ ρ ( (cid:126)r ) , n αβ , V KS ( (cid:126)r ) , λ αβ ] = − Tr ln iω + µ + ∇ − V KS − (cid:88) αβ ∈ L χ ∗ α ( (cid:126)r ) λ αβ χ β ( (cid:126)r (cid:48) ) − (cid:90) V KS ( (cid:126)r ) ρ ( (cid:126)r ) d r − λ αβ n αβ + (cid:90) V cryst ( (cid:126)r ) ρ ( (cid:126)r ) d r + E H [ ρ ( (cid:126)r )] + E LDAxc [ ρ ( (cid:126)r )] + Φ HF [ n αβ ] − Φ DC [ n αβ ] , (19)where n αβ is the occupancy matrix. Φ HF in Eq. (19) isthe Hartree-Fock approximationΦ HF [ n αβ ] = 12 (cid:88) αβγδ ∈ L ( U αγδβ − U αγβδ ) n αβ n γδ , (20)where indexes α, β, γ, δ refer to the fixed angular momen-tum L of correlated orbitals, and the matrix U αβγδ is theon-site Coulomb interaction matrix element. Similarly to DMFT, the on-site Coulomb interactionis already considered within LDA approximately, so it issubtracted, hence the double-counting term Φ DC . Oneof the popular choices is so-called “fully-localized limit”(FLL) whose form is Φ F LLDC [ n αβ ] = 12 ¯ U ¯ n (¯ n − −
12 ¯ J (cid:2) ¯ n ↑ (cid:0) ¯ n ↑ − (cid:1) + ¯ n ↓ (cid:0) ¯ n ↓ − (cid:1)(cid:3) , where ¯ n σ = (cid:88) α ∈ L n σαα , ¯ n = ¯ n ↑ + ¯ n ↓ , ¯ U = 1(2 L + 1) (cid:88) αβ ∈ L U αββα , ¯ J = ¯ U − L (2 L + 1) (cid:88) αβ ∈ L ( U αββα − U αβαβ ) . The constant matrix λ αβ is determined by the saddlepoint equations δ Γ LDA + U δn αβ = 0: λ αβ = δ Φ HF δn αβ − δ Φ DC δn αβ = (cid:88) γδ ( U αγδβ − U αγβδ ) n γδ − δ αβ (cid:20) ¯ U (cid:18) ¯ n − (cid:19) − ¯ J (cid:18) ¯ n σ − (cid:19)(cid:21) . (21)The FLL double-counting term tends to work quitewell for strongly correlated materials with very local-ized orbitals. However, for weakly correlated materials,FLL scheme describes the excessive stabilization of oc-cupied states and leads to quite unphysical results suchas the enhancement of the Stoner factor . In order toresolve the problems, “around mean-field” (AMF) wasintroduced in Ref. 67 and further developed in Ref. 66.One can say that the LDA+U method works whencorrelations are static - and at the same time local. Forexample cases where magnetic or orbital order are verywell developed. For a review of the LDA+U method seeRef. 68. Slave-Boson Method.
The physics of strongly cor-related electron materials requires to take into account -on the same footing, localized - quasi-atomic degrees offreedom, which are important at high energies, togetherwith strongly renormalized itinerant quasiparticles whichemerge at low energy. DMFT captures this physics via asophisticated quantum embedding that requires the solu-tion of a full Anderson impurity model in a self consistentbath. A less accurate but computational faster methdto solve the strong correlation problem, which precedesDMFT is the Gutzwiller method which has been shownto be equivalent to the slave boson method in the sad-dle point approximation. This approach starts with an1exact quantum many body problem, but one enlargesthe Hilbert space so as to introduce explicitly operatorswhich describe the different atomic multiplet configura-tions, and additional fermionic degrees of freedom whichwill be related to the emergent low energy quasiparticles.The method proceeded by writing a functional integralin the enlarged Hilbert space, supplemented by Lagrangemultipliers which enforce multiple constraints. The ap-proach, at zero temperature is very closely connected tothe Gutzwiller method, which appears as a saddle pointsolution in the functional integral formalism . In itsoriginal formulation this method was not manifestly ro-tationally invariant, but it was extended in this respect inRefs. 70–72. Further generalizations in the multi-orbitalformulation and to capture non-local self-energies was in-troduced in Ref. 73, and we denote this formulation asthe RISB (rotationally invariant slave boson) method.Within the slave particle method it is possible to go be-yond mean field theory, and fluctuations around the sad-dle point generate the Hubbard bands in the one particlespectra . The RISB method can be used to computethe energy of lattice models. When in conjunction withthe DMFT self-consistency condition, it gives the sameresults as the direct application of the method to the lat-tice model . In this review, we will restrict ourselvesto the RISB mean field theory, specifically from the per-spective of the free-energy functionals that describe thefree-energy of the system. We explain the physical mean-ing of the variables used in this method, and summarizesuccinctly the content of the mean field slave boson equa-tions using a functional approach. A precise operationalformulation of the method was only given recently . Forpedagogical reasons we start again with a Hubbard modelwith a tight binding one body Hamiltonian H ( (cid:126)k ).The variables used in RISB can be motivated by notic-ing that the many-body local density matrix ρ (with ma-trix elements (cid:104) Γ (cid:48) | ρ | Γ (cid:105) ) admits a Schmidt decomposition,which can be written in terms of the expectation valueof matrices of slave-boson operators φ Bn and φ † Bn , whichbecome φ Bn, φ ∗ Bn when the single-particle index α is M-dimensional, and can be stored as 2 M × M matrices Φ , [Φ] An ≡ φ An , [Φ † ] nA = φ ∗ An , so that: ρ = ΦΦ † . (22)The method also introduces Fermionic operators f α at each site (site indices are omitted in the following)which will represent the low energy quasiparticles at themean field level . The physical electron operator d isrepresented in the enlarged Hilbert space by d α = R αβ [ φ ] f β (23)where the matrix R , with elements R αβ at the mean fieldlevel, has the interpretation of the quasiparticle residue, relating the physical electron to the quasiparticles. When R is small it exhibits the strong renormalizations inducedby the electronic correlations. An important feature ofthe rotational invariant formalism is that the basis thatdiagonalizes the quasiparticles represented by the opera-tors f is not necessarily the same basis as that whichwould diagonalize the one electron density matrix ex-pressed in terms of the operators d and d † . Of centralimportance is the expression of the matrix R , in termsof the bosonic amplitudes: R αβ = (cid:88) γ (cid:88) ABnm F † α,A,B F † γ,nm [Φ † ] nA [Φ] Bm (cid:104) (∆ p (1 − ∆ p )) − / (cid:105) γβ = (cid:88) γ Tr (cid:2) Φ † F † α Φ F γ (cid:3) (cid:104) (∆ p (1 − ∆ p )) − / (cid:105) γβ . We introduced here the matrices F ,[ F α ] nm ≡ f (cid:104) n | f α | m (cid:105) f . The matrices∆ pαβ ≡ (cid:88) Anms (cid:104) m | f † α | s (cid:105)(cid:104) s | f β | n (cid:105) Φ † nA Φ Am = Tr (cid:2) F † α F β Φ † Φ (cid:3) have the physical interpretation of a quasiparticle densitymatrix: < f † α f β > = ∆ pα,β . For a multi-band Hubbard model with a tight-bindingone-body Hamiltonian H ( (cid:126)k ) and interactions Σ i H loci , theRISB functional, whose extremization gives the slave-boson mean field equations, is expressed in terms of φ i , φ † i (the slave-boson amplitude matrices) and the ma-trices λ ci , λ i , D . These are N × N matrices of Lagrangemultipliers: (i) λ ci enforces the definition of ∆ pi in termsof the RISB amplitudes (ii) λ i enforces the Gutzwillerconstraints and (iii) D i enforces the definition of R i interms of slave boson amplitudes. Another variable, E c enforces the normalization Tr[ΦΦ † ] = 1.The variables R , λ can be thought as a parametriza-tion of the self-energy. While the matrices D , λ c are aparametrization of a small impurity model (the dimen-sion of the bath Hilbert space is the same as that of theimpurity Hilbert space) , D is the hybridization functionof the associated impurity model while λ c parametrizesthe energy of the bath. ∆ p describes the quasiparticleoccupancies, which are the static analogs of the impurityquasiparticle Green’s function.The RISB (Gutzwiller) functional for a model Hamil-tonian with a local part which is bundled together witha local interaction term in H loc and a kinetic energy ma-trix which is non-local H ( (cid:126)k ) nonloc , was constructed inRef. 75:2Γ model [ φ, E c ; R , R † , λ ; D , D † , λ c ; ∆ p ] = − lim T → TN (cid:88) k (cid:88) m ∈ Z Tr ln (cid:32) i (2 m + 1) π T − R H ( (cid:126)k ) nonloc R † − λ + µ (cid:33) e i (2 m +1) π T + (24)+ (cid:88) i Tr (cid:20) φ i φ † i H loci + (cid:88) aα (cid:16) [ D i ] aα φ † i F † iα φ i F ia + H.c. (cid:17) + (cid:88) ab [ λ ci ] ab φ † i φ i F † ia F ib (cid:21) + (cid:88) i E ci (cid:16) − Tr (cid:2) φ † i φ i (cid:3)(cid:17) − (cid:88) i (cid:20) (cid:88) ab (cid:0) [ λ i ] ab + [ λ ci ] ab (cid:1) [∆ pi ] ab + (cid:88) caα (cid:16) [ D i ] aα [ R i ] cα (cid:2) ∆ pi (1 − ∆ pi ) (cid:3) cα + c.c. (cid:17) (cid:21) . (25)This method can also be turned into an ab-initio DFT+G method (or DFT+RISB). To motivate the con-struction of a
DFT+G functional we simply follow thesame path used above to go from the model DMFTHamiltonian to a DFT+DMFT functional. We substi-tute H ( (cid:126)k ) for − ∇ + V KS ( (cid:126)r ), which has a local anda nonlocal part, and follow the same steps as in theDFT+DMFT section.Γ DF T + G (cid:2) ρ ( (cid:126)r ) , V KS ( (cid:126)r ) , φ, E c ; R , R † , λ ; D , D † , λ c ; ∆ pi (cid:3) =Γ model [ H ( (cid:126)k ) → − ∇ + V KS ( (cid:126)r )] + Γ [ V KS ( (cid:126)r ) , ρ ( (cid:126)r )] − (cid:88) i Φ DC [∆ pi ]where Γ and Φ DC are the same functionals definedin the subsection on DMFT. The LDA+RISB and theLDA+G method are completely equivalent (more pre-cisely, the slave boson method has a gauge symmetry, anda specific gauge needs to be chosen to correspond to themulti-orbital Gutzwiller method introduced in Ref. 77.DFT+G was formulated in Refs. 78 and 79. The slaveboson method in combination with DFT was first usedin Ref. 80 in a non-rotationally-invariant framework andwith full rotational invariance in Ref. 73. For a recentreview see Ref. 81. Comparing the methods, critical discussion, future directionsand outlook
For weakly correlated systems we argued in section IIthat once the structure is known, we have a well-definedpath to compute their properties using DFT and the G W method. To go beyond requires to move in thespace illustrated in Fig 3. This has to be done whilerespecting as many general properties such as conserva-tion laws (Refs.82–84), sum rules, unitarity and causality(Refs.85–87) as possible. This is a very difficult problemwhich is under intensive investigation.This section reviewed several Green’s-function-basedapproaches available for studying strongly-correlated-electron materials. The reader may wonder why we con-sidered multiple methods. There are two reasons. First,as stressed throughout the paper and demonstrated in the examples, presented in the next sections there arematerials where correlations are mostly static, and oth-ers where dynamical correlations dominate the physics.These different types of correlations require differentmethods. Second, even when two methods treat the sametype of correlations, they have different accuracies andcomputational speeds. Finding the correct trade-off be-tween speed and accuracy will be important, in particularwhen high throughput studies start becoming feasible forstrongly correlated systems.As we strive towards a fully controlled but practicalsolution of the full many-body problem for solid statephysics, we will need more exact and thus slower meth-ods to benchmark the faster but more practical ones.Hence it is important to compare them and understandtheir connection. Static correlations can be treated byGW methods, and one can view the hybrid-functionalexchange-correlation potentials as faster approximationsto the QPGW exchange / correlation potential. One canalso assess whether the GGA (or LDA) exchange / cor-relation potential is a good approximation to the self-energy in a given material - by checking how close it is tothe corresponding self-consistent QPGW exchange cor-relation potential.In the same spirit one can understand the successes ofLDA+DMFT from the GW+DMFT perspective. Oneissue is the definition of U in a solid. The functional Φcan be viewed as the functional of an Anderson impuritymodel which contains a frequency-dependent interaction U ( ω ) obeying the self-consistency condition: U − = W − loc + Π loc . (26)This provides a link between LDA+DMFT, which usesa parameter U , and the GW+DMFT method, where thisquantity is self-consistently determined. An importantquestion is thus under which circumstances one can ap-proximate the Hubbard U by its static value. For projec-tors constructed on a very broad window, U ( ω ) is con-stant on a broad range of frequencies . An importantopen question is how one can incorporate efficiently theeffects of the residual frequency dependence of this inter-action.Another question is the validity of the local ansatz forgraphs beyond the GW approximation. This questionwas first addressed in Ref. 89, who showed that the low-est order GW graph is highly non-local in all semicon-3ductors, which can be understood as the exchange Fockgraph is very non-local. On the other hand, higher-ordergraphs in transition metals in an LMTO basis set wereshown to be essentially local.Consider a system such as Cerium, containing light spd electrons and heavier, more correlated, f electrons. Weknow that for very extended systems, the GW quasipar-ticle band structure is a good approximation to the LDAband structure. Therefore the self-energy of a diagram-matic treatment of the light electrons can be approxi-mated by the exchange-correlation potential of the LDA(or by other improved static approximations if more ad-mixture of exchange is needed) . Diagrams of all ordersbut in a local approximation are used for the f electrons.In the full many-body treatment Σ ff is computed us-ing skeleton graphs with G loc and W loc . To reach theLDA+DMFT equations, one envisions that at low en-ergies the effects of the frequency dependent interaction U ( ω ) can be taken into account by a static U , whichshould be close to (but slightly larger than) U ( ω = 0).The f f block of the Green’s function now approachesΣ ff − E dc .We reach the LDA+DMFT equations with some addi-tional understanding on the origin of the approximationsused to derive them from the GW+DMFT approxima-tion, as summarized schematically inΣ GW + DMF T (cid:16) (cid:126)k, ω (cid:17) −→ (cid:32) ff − E dc (cid:33) + (cid:32) V xc [ (cid:126)k ] spd,spd V xc [ (cid:126)k ] spd,f V xc [ (cid:126)k ] f,spd V xc [ (cid:126)k ] f,f (cid:33) . Realistic implementations of combinations of GWand DMFT have not yet reached the maturity ofLDA+DMFT implementations, and are a subject of cur-rent research. Recent self-consistent implementations in-clude Refs. 60 and 90.When strong dynamical correlations are involved, thespectra is very far from that of free fermions. The oneelectron spectral function A ( (cid:126)k, ω ) displays not only a dis-persive quasiparticle peak, but also other features com-monly denoted as satellites. The collective excitationspectra, which appear in the spin and charge excitationspectra, does not resemble the particle-hole continuum ofthe free Fermi gas with additional collective (zero sound,spin waves) modes, produced by the residual interactionsamong them. Finally the damping of the elementary exci-tations in many regimes does not resemble that of a Fermiliquid. Strong dynamical correlations are accompaniedby anomalous transport properties, large transfer of op-tical spectral weight, large mass renormalizations, as wellas metal-insulator transitions as a function of tempera-ture or pressure. These can be captured by DMFT, whichcombined with electronic structure, enable the treatmentof these effects in a material-specific setting, but not byLDA+G which only provides a quasiparticle descriptionof the spectra. Many successful comparisons with ex-perimental ARPES and optical and neutron scattering data have been made over the last two decades usingLDA+DMFT which makes an excellent compromise ofaccuracy for speed, and it is now the mainstay for the elu-cidation of structure property relations in strongly cor-related materials. LDA+G can only describe at best thequasiparticle featurs in that spectra.On the other hand, as it will be stressed through ex-amples, for total energy evaluations - which are a cen-tral part of the material design workflow, faster methodsare currently needed. We described above two methods,the LDA+U method, and the Gutzwiller RISB method,which fall in this “fast but less accurate” category. Thesemethods can be viewed as approximations to the manybody problem within a DMFT perspective. As pointedout in Refs. 75 and 91, the Gutzwiller RISB leads toa DMFT-like impurity solver with a bath consisting ofonly one site. LDA +U can be viewed as a limiting caseof DMFT, where a static local self-energy is considered.There are numerous algorithmic challenges in optimiz-ing studies of materials based on DMFT. While CTQMCruns for solving the Anderson impurity model, i.e. thesingle orbital case, as well as 3 or 2 orbitals ( t g and e g electrons) can be completed on one CPU in less thanone day for extremely low temperatures, a full d-shell(5 electrons) requires several days, and the full f-shell isstill at the border of what can be done with current meth-ods. All this assumes high symmetry situations, wherethe hybridization function is diagonal. Off-diagonal hy-bridization introduces severe minus sign problems. Al-ternative exact diagonalization-based methods, such asNRG or DMRG will be needed. This would also helpwith the problem of reducing the uncertainties involvedin the process of analytic continuation.While the ansatz 10 has reproduced the photoemis-sion spectra of many materials, there have not been high-throughput studies which would enable us to systemat-ically search for deviations. This requires the improve-ment of computational tools, an area of active research.What if the k and ω dependencies cannot be disentan-gled? This situation may arise near a quantum criticalpoint. Methods to incorporate the non-local correlationsbeyond DMFT are an important subject of active re-search, which is reviewed in Ref. 59.Armed with an understanding of methods to treat cor-relations and their physical and computational trade-offs,we proceed in section IV to construct a workflow for de-signing correlated materials. IV. MATERIALS DESIGN WORKFLOW
Condensed matter physics has a long standing tradi-tion of constant interplay between theory and experimen-tation. The field of strongly correlated electron systems,has been driven by unexpected experimental discovery anew materials, followed by a large number of theoreti-cal ideas which get refined as new experimental informa-tion becomes available. This is described in Fig. 9 panel4
Qualitative IdeasPhenomenologySimple model Calculation ExperimentationTheoryAlogrithmsMaterial-specificComputationQualitative IdeasPhenomenologySimple model Calculation Experimentation (a)(b)
FIG. 9. (a) Traditional theory experiment theory interaction.This loop is initiated by an experimental discovery for whichmany theories are proposed and corrected by further exper-imentation. (b) The modern material design loop augmentsthe theoretical contribution, now this loop can be initiated byeither theory or experiment. (a). Panel (b) describes how theory, algorithms and com-putational power have enhanced theoretical capabilities,which make the approach material-specific, thus enablingtheory-assisted material design.At the current state of development, the material de-sign process for strongly correlated electron compoundsshould begin at the intuitive level, for example somephysical idea of model that one would like to explore ortest, a property of a strongly correlated material whichcould result in a useful property, a class of compoundsone would like to investigate comparatively, some ideas ofchemistries of strongly correlated materials which wouldenhance desirable solid state properties. This zeroth or-der step can be refined with simplified quantitative cal-culations using model Hamiltonians, or other computa-tional tools which refines our intuition. After this zerothorder step, it is natural to divide the process of materi-als design into three additional steps as summarized inTable II, which will be illustrated by example in the fol-lowing sections of this paper.The first step is the quantitative calculation of the elec-tronic structure , namely, how to go from a well definedstructure (i.e. atomic positions and ionic charges) tophysical or chemical properties. Given a crystal struc-ture, we seek to compute one or several electronic prop-erties such as orbital occupancies, transport coefficients(resistivities, mobilities, thermoelectric coefficients, etc.)Mott or charge transfer gap sizes, magnetic order param-eters, etc. As we have seen in section II, for correlatedelectron materials the computational method used forelectronic structure depends on the strength and kind of correlations and should be chosen appropriately. Some-times, this step is divided into two. In the first step firstone derives from first principles an effective model Hamil-tonian, and then one solves the model and explores itsconsequences. While the advanced functionals describedin the previous section go directly from structure to prop-erty bypassing the model Hamiltonian, the latter can beuseful for the interpretation of the results.The second step is structure prediction : predict thecrystal structure given a fixed chemical composition. Asuccessful prediction would take a formula, like Fe O forexample, and return the correct crystal structure–in thiscase, the composition forms in the corundum structure,called the α phase. For a more complete characteriza-tion, we would seek to predict not only the ground statestructure, but nearby local minima as well, termed poly-morphs. Again taking Fe O as an example, this compo-sition also forms in the spinel structure, found naturallyas the mineral maghemite and termed the γ phase, aswell as cubic β and orthorhombic (cid:15) phases. Polymorphsgenerally are formed in different temperature and pres-sure regimes, and modeling these effects add an addi-tional layer of complexity. However, simply enumerat-ing the low-energy local minima at zero temperature canalready provide a broader picture of the structural di-versity of a composition. Furthermore, if the structureone is interested turns out not to be the ground statebut a metastable structure one can design methods forstabilizing it, either by choosing an appropriate synthe-sis method or by applying external perturbations such asstress exerted by a properly chosen substrate.The general procedure for structure prediction involvesplacing atoms in a unit cell and using an algorithm to ef-ficiently traverse the space of atomic configurations andcell geometries to arrive at low energy structures. Thisstep requires having an accurate method for producingthe energy of a given configuration of atoms and sam-pling these configurations. There are a number of struc-ture prediction techniques (see the review Ref. 92) in-cluding the simulated annealing approach , evolution-ary algorithm methods , structure models by anal-ogy based on data mining and machine learning ,metadynamics , basin and minima hopping ,random structure searching , and so on.The third step is testing for global stability : given thelowest energy structure of a fixed composition, checkwhether it is stable against decomposition to all othercompositions (phase separation) in the chemical system.The steps involving total energies are assisted by elec-tronic structure methods and material databases. In par-ticular the third step which requires the knowledge ofall other known stable compositions, their crystal struc-tures and total energies, is now facilitated by materialsdatabases containing data in standardized computableformats, such as the Materials Project , the Open Quan-tum Materials Database (OQMD) , AFLOWlib andNIMS . With this information, the energetic convex hullfor a chemical system can be constructed and the target5 Name Step Tools
Motivation and analysis defining directions and hypotheses to be tested,
Heuristics, theory, simplified computations refining them with simplified computations electronic structure structure → property electronic structure codes, DFT, DFT+DMFTstructure prediction composition → structure evolutionary algorithms, Monte Carlo, minima hoppingglobal stability chemical system → composition convex hull from materials databasesTABLE II. Three step workflow of materials design. Electronic structure and structure prediction have been accessible fordecades, with density functional theory (DFT) as the key underlying method. In contrast, checks for global stability requireknowledge of all other structures within a chemical system, which was only possible after the creation of extensive computationalmaterials databases. composition checked for stability against decomposition(phase separation).For weakly-correlated materials, the entire workflowcan be built around LDA/GGA for total energies and G W for spectral properties. For correlated materials,GGA is a good starting point for computing total en-ergy differences (such as reaction energies or structuralenergy differences). In this review we highlight some fail-ures of the LDA/GGA energies in structural predictionand phase stability, and ways to introduce corrections toaccount for the correlation effects.The need to correct DFT total energies for materialsdesign projects, is broadly recognized in the context ofall the material databases where the GGA/LDA resultsare corrected using semi-empirical schemes. There arethree broadly used schemes in literature, and they haveassociated databases, the fitted elemental-phase referenceenergy (FERE) scheme , Materials Project , and OpenQuantum Materials Database (OQMD) .While the details of the implementation are different,they have two key elements in common. First, insteadof using GGA, the total energies are computed withinthe GGA+U method, with some U values empiricallyassigned to each element. Second, the experimental for-mation energies ∆ H exp are used to determine best fits forelemental energies E Fitted ( A ), where A is an element, fora training set of compounds by solving the linear least-squares problem.∆ H exp (A m B n ) ≈ E GGA+Ucorrected (A m B n ) (27)= E GGA+U (A m B n ) − mE Fitted ( A ) − nE Fitted ( B )We note that all elemental energies are fitted for allelements in FERE (whithin the set of relevant elements)while only selected ones are fitted in MP and OQMD(especially in the “fit-partial” scheme in OQMD).We describe the different correction methods in Ap-pendix A. In this article we will use the Materials Projectdatabase for analysis of phase stability and estimateprobabilities based on their data. We will show in sectionVI E that careful consideration of correlations is essentialnot only for calculation of phase stability, but also forstructure prediction.Notice that the theoretical workflow, outlined in Ta-ble II progresses in an order different from experimen- tal solid state synthesis. There, elements and simplecompounds in a chemical system are combined and sub-jected to heating/cooling cycles to provide the kineticenergy necessary for atomic rearrangement to form newstoichiometries (of which there may be more than one).Finally the stoichiometries crystallize to form structureswhich are then isolated for the study of their properties. V. STATISTICAL INTERPRETATION
DFT has reached a high degree of stability and scal-ability, enabling software packages such as USPEX to implement genetic algorithms on top of DFT to suc-cessfully predict never before observed structures. As dis-cussed in the previous sections, correlations in the form ofU and empirical corrections are now available in severaldatabases. Since these methods are not exact methodsand suffer from systematic errors, compounds predictedto be stable will not necessarily be found in experiment,and vice versa.The main question we address in this section is theinterpretation of the above-hull/below-hull energies thatwe compute within LDA/GGA (with or without the em-pirical corrections). There are two related questions: (1)what is the likely error in the computed energy, (2) howlikely is the compound to be synthesized given its energyrelative to the convex hull. Namely, how likely is oneto find the target compound? This assessment serves asa background for the conclusion section VII, where weevaluate the results of various material design projects.Reference 107 modeled the computational error - thedifference between computed and experimental formationenergies - as a random variable with normal distribution.A normal-distribution was also used by Ref. 4, as seen inFig. 11. We follow a similar statistical approach to thequestion.Denote by E exp the experimental heat or enthalpy peratom of the reaction A+B → AB at low temperature,and denote by E calc the same quantity computed usingan approximate method, like GGA, or the empirically-corrected value of this computation. We treat E calc and E exp as real random variables, and analyze the distribu-tion of the variable d = E calc − E exp :6 P ( E calc − E exp = d ) = F αβµ ( d ) (28)where F αβµ ( d ) is some probability distribution func-tion (PDF) with center µ and scale α , as well as someshape parameter β . Since GGA (or one of its correctionschemes) is reasonably accurate, we expect F αβµ ( d ) tobe concentrated around the center.In order to study this distribution, we observe that thesame quantity in Eq. (28) describes computational errorsin energies above-the-hull as well as computational errorsin formation energies. This holds because the distribu-tion applies to computational error in reaction energiesin general . The crucial point that makes this possible isthat the number of atoms is balanced on the left and onthe right (so as to cancel core energies). Therefore wecan train a statistical model on the distribution of com-putational error for formation energies, for which thereexists a reasonably-sized experimental data set, and thenmake predictions based on above-hull energies.We expect a stronger statistical-correlation when allthe systems A, B, and AB are weakly correlated thanwhen the correlations are strong, hence the parameters α, β, µ should be taken in a well-defined space of materi-als defined from the outset. Since we will use this modelfor prediction, it is necessary to fit the parameters usinga large-enough sample of representative materials.Figure 10 shows the distribution of computational er-ror E calc − E exp for formation energies of compoundscollected in the data set. The data set includes 1,500substances from the OQMD database and the materialslisted in table II in [14]. It is noteworthy that the exper-imental formation energies are available as well as crys-tal structures of the compounds in the OQMD database(query was used to collect these pairs). The experi-mental data originates from 2 sources: the SGTE SolidSUBstance (SSUB) database , and the thermodynamicdatabase at the Thermal Processing Technology Centerat the Illinois Institute of Technology (IIT) , as de-scribed in Ref. 4. The experimental formation energieswere also collected from table II in 14 and were mergedwith the data set. In Fig. 10 the left-side plots showthe computational error for pure GGA calculations (re-produced with our own VASP calculations for the com-pounds in the data set), whereas the right hand side in-cludes +U for some of the compounds, as in the fromthe Materials Project’s recipe (see appendix for details).We observe that the distribution is skewed to the right,and that application of the U correction is not enough toundo the skewness, although it reduces σ from 427 meVto 266 meV (with an increase of the mean error from 155meV to 172 meV). This data is summarized in Table III.As discussed above, various authors have corrected theGGA/GGA+U values for formation energies by shift-ing the chemical potentials of elements. This is usuallyachieved by fitting the set of experimental results on theequations for GGA/GGA+U formation energies as in Eq.(27). The effect of this procedure on the distribution F αβµ ( d ) is to make it centered, and to eliminate the driftof the bivariate distribution. This can be seen in theOQMD fit (Fig. 11, as well as in FERE and MaterialsProject (Fig.15).Another important property of the corrected distribu-tions is that the error can be seen as approximately in-dependent of the value of the computed energy: P ( E calc − E exp = d | E calc = x ) ≈ P ( E calc − E exp = d )= F αβµ ( d )This (approximate) independence is evident in Fig. 12.One can reason that as long as the computational erroris small, it should not be correlated with the value ofthe computed energy. Finally, with just a few 1000’sof points, there is not enough data to split the domaininto sub-ranges and make meaningful statistical analysis.With more data one could refine the distribution param-eters on ranges of E calc = x , or possibly other variablesthat we did not consider here.For prediction, the probability that AB forms as theground state , when the computed energy is at a distance x above the Hull is given by P ( x ) = (cid:90) −∞ P ( E exp = y | E calc = x ) dy (29)= (cid:90) −∞ P ( E calc − E exp = x − y | E calc = x ) dy ≈ (cid:90) −∞ F αβµ ( x − y ) dy = 1 − F αβµ ( d (cid:54) x )where F αβµ (d (cid:54) x) is the cumulative distribution func-tion corresponding to F αβµ . This expression can beevaluated numerically by estimating the number of datapoints in the distribution tail (where P ( x ) ≈ { ( E calc , E exp ) | E calc > x } { ( E calc , E exp ) } , however, it is convenient to have an analytic form. Forthe corrected distributions, we postulate that F αβµ isNormal or generalized-Normal distribution - which in-cludes also Normal ( β = 2), exponential ( β = 1), as wellas uniform ( β = ∞ ) distributions: F αβµ ( d ) = β α Γ(1 /β ) e − ( | d − µ | /α ) β (30)Using the experimental data (see above), we calculatedthe numbers listed in Table III. We used the maximumlikelihood method to estimate α, β and µ . The first col-umn summarizes the distribution of raw GGA calcula-tions, which were reproduced with our own GGA runs forthe compounds in the data set. The second column cor-responds to Materials Project’s GGA or GGA+U (GGA7 FIG. 10.
Left : Distribution of computational error reproduced with VASP GGA(PBE) for compounds collected in the data set.The data set includes 1,500 substances from OQMD and the materials listed in table II in 14. To get the experimental formationenergies for the corresponding compounds in the data set, we use the OQMD database, where the experimental data is availableas well as the crystal structure. The experimental data is also merged with table II in 14. The pure GGA(PBE) formationenergies are calculated for the compounds in the data set and compared with the corresponding experimental formation energies.
Right: plots for Materials Project’s raw data, which adds nonzero U for some elements in certain compounds. As can beseen, addition U is not enough to make the distribution F αβµ ( d ) un-skewed or fix the drift to the bottom (this is also evidentin the Fit-none distribution from OQMD , shown in Fig. 11). Also shown are the ellipses of the bivariate-normal estimators,demonstrating the drift. for most compounds and GGA+U only for certain cor-related materials specified in appendix), which was eval-uated using the raw data from Materials Project. Thesedistributions are depicted in Fig. 10. We did not includemodel estimates for the non-corrected distributions, sincethey do not comply with some of the assumptions (as ex-plained above). The third column in the table is fittedfor ∆ H F ERE from Table II in Ref. 14. As expected,the parameters for FERE show a smaller standard er-ror (for a smaller dataset). The fourth column corre- sponds to corrected formation-energies from the Materi-als Project. Again, the parameters show a smaller stan-dard error compared to bare GGA. These distributionsare depicted in Fig. 15. Figure 14 shows the calculatedprobabilities in each one of the schemes for E calc .Interestingly, the mean values for pure GGA and Mate-rials Project’s raw data (table III) are positive, meaningthat GGA/GGA+U tend to over-estimate formation en-ergies. We also observe that the corrected distributionsare more exponential-like then normal. In fact, they all8 FIG. 11. Comparison of “Fit-None” to “Fit-Partial” inOQMD (from Ref. 4). Left side shows raw data (with Uadded to some correlated elements in certain compounds) -there is an evident skewness and drift. On the right the datais corrected by a small set of chemical potentials. The skew-ness and drift are eliminated.FIG. 12. Plot of the computational error in Materials Projectformation energies against the calculated formation energy,demonstrating little or no correlation between E exp −E calc and E calc . Each data point corresponds to a pair of calculated andexperimental formation energy for a material in the data set. get a very low score on statistical Shapiro testing forNormalcy, but pass the Kolmogorov test for generalized-Normal distribution.In general we would like to apply the same kind of sta-tistical error analysis to the systematic computational er-ror made by correlated electronic structure calculations,such as DMFT, or Gutzwiller. However, since there ex-ists no database of energies for these methods, at thispoint in time it is impossible to estimate the parame-ters of the model Eq. (30) for correlated methods. Still, FIG. 13. The determinant reaction for a new material X- X lies either above or below the convex hull of reactions.In the illustration above it is well below the ABC triangle,although its energy is above the B compound (note that inthis case X still forms out of A, B, and C). X’s computeddistance below the hull, which we denoted in the text by E X , isthe quantity that determines its formation probability in ouranalysis. In this diagram it is equal to X’s vertical distancefrom the triangle ABC, and therefore the reaction ABC → Xis the determinant reaction in this hull.FIG. 14. P ( E calc ), probability for a compound to exist, giventhe computed formation energy E calc with the two correctionschemes summarized in Table III. since we expect these methods to be more accurate thanGGA, we expect their standard error σ to be smaller,and their mean µ to be closer to 0 than what we foundfor GGA. Since our probability estimate depends onlyon the distance from the convex hull, it is determinedonly by one triangle in the diagram, which we refer toas the determinant reaction, see Fig. 13. Only the en-ergies of the relevant compounds (in the diagram, theseare A, B, C, and X) need to be re-estimated using thecomputationally expensive correlated method, and thusone can improve the probability estimate for X, withoutrequiring too much extra computation. Since the cumu-lative distribution function (Eq. (29)) for the correlatedcomputation will be sharper, if the energy of X falls be-low the hull, its probability will clearly be larger than0.5, whereas if the energy is well above the hull we arguethat it is unlikely to exist.9 FIG. 15. Distribution of computational error for empirical correction schemes: FERE (left) and Materials Project correctedscheme (right). For FERE, the data comes from table II in Ref. 14. For the Material Project, the corresponding corrected GGAvalues were downloaded using their Python interface. Fitting of the distribution was done using the maximum likelihood method fit() (from Scipy ) for the generalized-normal distribution ( scipy.stats.gennorm ). The distributions are summarized in thethird and fourth columns of table III.
Narayan et al. constructed a related statisticalmodel to address theoretical predictions of new materi-als using formation energies computed with the MaterialsProject corrections. They proposed an optimal cutoff en-ergy (cid:15) = 0 . eV , to minimize the rates of false-positivesand false-negatives. They suggested that the synthesis ofmaterials with (cid:15) < (cid:15) should be attempted. Our analysiscan be seen as a refinement of the idea of a cutoff.From an experimental standpoint, the question to an-swer is whether an experimental search should be pursuedor not. This depends not only on the function P ( x ) butalso on how interesting is the estimated property that weare seeking. Also, notice that we are asking the questionof whether we will obtain a material close to its ground state. It is well known that many known compounds aremetastable . Therefore, when looking for compoundswhich have not been synthesized before, we are obtaininga lower bound for the probability of finding a new ma-terial, which is in greater than P ( x ). Third, notice thatstability analyses based on convex hulls can only ruleout the stability of a compound, rather than definitivelyprove its stability, since technically there are an infinitenumber of compositions and structures to consider in thephase diagram.Theory assisted material discovery has already hadnotable successes, where the predicted compounds weresuccessfully synthesized. It has lead to new multiferroicmaterials (Refs. 7 and 114). The material design strat-0 Before corrections After correctionGGA (U=0) MP (+U) FERE MPMean 136 151 6 -31MAE 155 172 51 130 σ
427 266 81 192center µ – – 1 -16shape β – – 0.85 0.90scale α – – 39 109 F . Maxi-mum likelihood values for the parameters of the distributionfunction in Eq. (30) are given for the correction schemes.All units are meV except for β (unitless) and the number ofpoints. MP (+U) stands for Materials Project’s raw data,which only includes the U correction for certain correlatedcompounds, as described in appendix. egy was used to discover and synthesize 18 new ternarysemiconductors . It has also been used to find two newstructures in the Ce-Ir-In system . There are also reportsof broad theoretical searches where experimental synthe-sis did not find the theoretical approach predictive .As an example from that work, we examine KScS whichhas E calc = − . eV relative to the convex hull. Ac-cording to our analysis it has probability 0.8 to exist andindeed it has been synthesized successfully in Ref. 115. VI. CASE STUDIES
In the following, we illustrate the concepts of the previ-ous section providing examples of the design of correlatedmaterials. We limited ourselves to a few projects we arefamiliar with and omitted important areas, as for ex-ample the design of nuclear fuels , where strongcorrelations in the solid state play a major role.The examples encompass materials displaying super-conductivity VI A, VI B, VI D, charge disproportionationVI C, and metal-insulator transitions VI E.The motivation for material design projects in the fieldof strongly correlated electron materials is not only tofind new compounds , but to use the process to check ourunderstanding of the theory, validating correct ideas andapproaches and discarding wrong ones. In each case thematerial design project begins with some intuitive ideathat one would like to test. We then move through theprocedure outlined in section IV , describing the steps togo from structure to property, from composition to struc-ture and from components to composition. We concludeeach section with the lessons learned from the project.We highlight the tools currently in use and the role thatcorrelations play in each stage of the workflow.The field of material design of strongly correlated elec-tron systems is in its infancy and the examples spanthe development of the workflow itself—in some cases inthe examples, only part of the full workflow was applied in the design process. Throughout the presentation westress the importance of qualitative ideas and chemicalprinciples and how they can be supported and enhancedwith modern computational techniques.
A. Tuning the charge-transfer energy
Background
The cuprate superconductors are classic examples ofcorrelated materials. Since this family of compounds ex-hibits the highest known superconducting transition tem-peratures at ambient pressure (surpassed only by H S athigh pressures), they have been intensely studied sincetheir discovery in the mid-1980s. They have motivatedan immense body of work in condensed matter physics,both in new theories and huge leaps in experiment.Structurally, all the cuprate families have in commonCuO planes which support superconductivity. They aredescribed by the chemical formula XS n − (CuO ) n , where n CuO planes are interleaved with n − T c ( T c,max )occurs. Empirically, it is known that T c,max is stronglymaterials-dependent, ranging from 40 K in La CuO to138 K in HgBa Ca Cu O .It is now agreed that superconductivity arises fromdoping a parent compound which is a charge transferantiferromagnetic insulator , and that the symmetry ofthe superconducting state in these materials is d -wave,as it was predicted theoretically . One shouldstress however, that there is no consensus on the mech-anism that control the superconducting critical temper-ature and as a consequence what leads to the correla-tion of T c,max with the apical oxygen distance. Even ifone accepts a mechanism, such as the proximity to mag-netism and the presence of on-site Hubbard correlations,controversies persist, as both weak coupling and strongcoupling approaches predict the correct nature of the su-perconductivity. Setting the Targets & Framing the Questions – Couldmaterials design help us verify or falsify theories of thehigh-temperature cuprate superconductors? Or at leastnarrow down the set of competing theories and sharpenour theoretical understanding while suggesting new com-pounds in this important class of materials?We chose the T -type layered perovskite La CuO asthe starting point. Our intuition led us to propose thesite substitution of the apical oxygen with sulfur. Due tothe larger ionic radius of sulfur as compared to oxygenwe expect the LaS charge reservoir layer to be crowded.To compensate, we explored the effect of substitutingthe large La ion with smaller trivalent ions R , selectedfrom the lanthanide-like elements. The compositions1 FIG. 16. Cluster DMFT investigation of a model ofLa − x Sr x CuO (Weber et al. ). Zero temperature phasediagram connects structure to physical properties (supercon-ductivity and antiferromagnetism). we considered were R CuS O (shown in Fig. 19) and R CuSO . We include the monosulfide in hopes that theconfigurational entropy of only replacing a quarter of theapical oxygens with sulfur would help stabilize the targetphase. Structure to Property
Cluster DMFT can describe the competition and syn-ergy between antiferromagnetism and d -wave supercon-ductivity .Weber et al. used Cluster LDA+DMFT with an ex-act diagonalization to bridge between the structure ofLa − x Sr x CuO (LSCO) and the physical properties ofmagnetism and superconductivity. As shown in Fig. 16,the method describes well the zero temperature phase di-agram of this compound, including the antiferromagneticphase (characterized by the magnetic moment S z ) andthe d -wave superconductor characterized by the anoma-lous self energy and the superconducting gap ∆.Multiple studies concluded that T c,max is an increasingfunction of the apical oxygen distance . Changingthe apical oxygen position, however, affects many differ-ent parameters of the electronic structure of the copperoxygen planes, ranging from the hoppings to the chargetransfer energy and to the relative positions of the d x − y and d z orbitals. For a working design process, we need toelucidate the connection between these parameters and T c,max .To investigate the link between T c,max and the variousparameters describing the electronic structure, one needsto vary them independently, which is easy to do in-silico,based on a combination of first principles calculationsand dynamical mean-field theory to directly model thesuperconducting state. Specifically, La − x Sr x CuO hasthe largest charge transfer energy and a small T c,max , we vary the charge transfer energy and the hoppings inde-pendently and the results are shown in Fig. 18. It isclear that (within LDA+DMFT) the superconductivityincreases rapidly when we decrease the charge transfergap, while it is very weakly sensitive (and in fact de-creases) when we increase the oxygen oxygen overlap andthe hopping integral t (cid:48) (depicted in Fig. 17) on the squarelattice (not shown).Inspired by these results and the earlier heuristic con-siderations of connecting structure to property, we aimedto design new cuprates with reduced charge transfer gaps(and thus higher T c ) via sulfur (S) substitution. Tocheck that S substitution helps control the charge transfergap, and to see which other variables it affects, we usedthe maximally localized Wannier function imple-mented in the WANNIER90 code .The results for the GGA estimates of the charge trans-fer energies are displayed in Fig. 19 and the secondcolumn of Table V together with the reference systemLa CuO . Both figure and table suggest that Sc CuS O having the highest charge transfer energy would have thesmallest T c,max among proposed cuprates R CuS O andLa CuS O having smaller charge transfer energy thanLa CuO would have higher T c,max than that of La CuO - based on the charge-transfer theory in Ref. 126.Sakakibara et al. observed a positive correlationbetween T c,max and the crystal field splitting ∆ E = (cid:15) d ( x − y ) − (cid:15) dz based on a two-band model incorporatingboth d x − y and d z orbitals. They found that large crys-tal field splitting ∆ E , which is also related to the reduc-tion of d z contribution to the Fermi surface, enhances T c,max . In order to test their theory, we examine ∆ E inFig. 20 and the last column of Table V. The crystal fieldsplitting ∆ E in the oxysulfides family R CuS O is gen-erally smaller than that in La CuO , implying that sulfursubstitution would suppress T c,max . La CuS O havingthe smallest crystal field splitting ∆ E among the oxysul-fides family, for example, has the smallest T c,max basedon the theory, while the charge-transfer theory gives theopposite result. Therefore a detailed study of La CuS O is necessary to test the validity of these two theories. Structure Prediction – To check for local stability weused a 2 × × T (cid:48) -type layered perovskite, know-ing that substitution of the large La ion for the smallerPr and Nd led to a rearrangement of the charge reser-voir layer into the fluorite structure. We found that the T -type structure was indeed stable and there was no out-of-plane buckling, although the CuO octahedra favoredaxial rotations ( a a c − p in Glazer notation). See Fig. 21. Global stability – We checked the thermodynamic sta-bility of the proposed compounds against competingphases by selecting commonly known reactants and com-puting the formation enthalpies of the synthesis path-ways as shown in Table IV. We computed the total en-ergies of formation ∆ E = E products − E reactants , andfind that all differentials are positive, indicating the re-2 HgBa CuO Hg(CaS) CuO FIG. 17. The empirical correlation of T c,max with the apical oxygen can be due to variation of many parameters in the lowenergy Hamiltonian describing different processes (top left), t (cid:48) /t as originally suggested by Pavarini et al. (top right), as wellas the charge transfer gap as suggested by Weber et al. (bottom) . We have updated the bottom figure to include Hg-basedcuprates.FIG. 18. Cluster DMFT investigation of a model of the copper oxide planes (Weber et al. ). After showing that the modeldisplays the correct phase diagram (top), starting with La − x Sr x CuO (LSCO) one varies independently the parameters thatcontrol the charge transfer energy (cid:15) p − (cid:15) d (bottom left) and the oxygen oxygen overlap t pp (cid:48) (bottom right). The figure displaysthe dependence of the maximum of the superconducting order parameter (which serves as a proxy for T c,max ). We see that themaximum superconducting gap ∆ max is very sensitive to the charge transfer gap, which thus controls (in this strong couplingtheory) the critical temperature and not with t (cid:48) or t (cid:48) pp . ε d – ε p ( e V ) d apical (Å) Ga CuS O La CuO Sc CuS O Lu CuS O Y CuS O La CuS O HgBa CuO Hg(CaS) CuO FIG. 19. Explored compositions, generated by substitutingthe apical oxygens (arrows) and rare earth ion (green spheres)in La CuO to form the family of compounds R CuS O .Two mercury compounds HgBa CuO and Hg(CaS) CuO are also added for comparison. Plotted are the extractedcharge-transfer energies vs. the apical distance of sulfuratom in the proposed compounds after structural relaxationin LDA. Note that charge-transfer energies in known cupratessuperconductors are observed in the range from ∼ . ∼ . actions target phases are unfavorable. However, it isknown that many functional materials are metastable,protected from decay by large energetic barriers. Theparent cuprate La CuO is an example: as shown on thelast line of Table IV, La CuO is actually unstable by65 meV/atom. We also examined the volume differen-tials ∆ V = V products − V reactants with the knowledge thatoften high pressure synthesis allows otherwise unstablecompounds to form. Notice that ∆ V are overwhelmingnegative, meaning the application of high pressure mayallow the formation of the target phases.To look at the question of global stability ofLa CuS O and La CuSO , we use the modern mate-rials databases, and reanalyze the entire La-Cu-S-O sys-tem to construct the convex hull (plotted in Fig. 22) andglobally investigate stability.Using the convex hull, we can assess the stability ofthe reactants and products reported in experiment. InFig. 23, we plot the energies relative to the convex hullfor all reported compounds. Negative values are stabilityenergies against decomposition. We find that La CuS O and La CuSO are unstable at 232 and 324 meV/atomabove the hull, respectively, which is slightly more un-stable than the earlier estimates which only tested a fewreactions, and we gain additional information as we learnthese compounds in equilibrium would decompose into:La CuS O → La SO + CuS4La CuSO → SO + 4Cu + La SO Additionally, LaCuSO lies 17 meV/atom below the ∆ E - b a nd ( e V ) d apical (Å) La CuO Sc CuS O Lu CuS O Y CuS O La CuS O Hg(CaS) CuO HgBa CuO FIG. 20. The crystal field splitting (cid:52) E − band = (cid:15) d ( x − y ) − (cid:15) dz computed by downfolding into the effective two-bandmodel. La CuO and the copper oxysulfide family areshwon with two mercury compounds HgBa CuO andHg(CaS) CuO for comparison. (cid:52) E − band in the oxysul-fides is generally smaller, implying stronger mixing between d z and d x − y orbitals. Thus sulfur substitution would sup-press T c,max according to the two-band theory , while T c,max would be enhanced according to the charge-transfer theory .Adapted from Ref. 134.FIG. 21. Octahedral rotations in Sc CuS O as shown by asection of the CuO plane with four atoms. Adapted fromRef. 134. hull, so it is likely to be stable according to the anal-ysis of section V.Following our proposal in section V, we obtainsharper probability estimates using a correlated method- GGA(PBE) + Gutzwiller. As outlined, we only need toexamine the energies of the determinant reaction. Theon-site Coulomb interaction U = 8 eV and Hund’s cou-pling constant J = 0.88 eV are used in the Gutzwillercalculations . We find that La CuS O and La CuSO are still unstable at 214 and 286 meV/atom abovethe hull, respectively - consistent with a recent experi-ment . Conclusion – Experimental support for the idea thatthe reduction of the charge transfer gap results in an4 ∆ E ∆ V Synthesis pathway232 -1.92 La SO + CuS → La CuS O
362 -1.45 Y SO + CuS → Y CuS O
415 -1.36 Lu SO + CuS → Lu CuS O
542 -1.04 Sc SO + CuS → Sc CuS O
304 -2.00 La O + CuS → La CuSO
780 -1.11 Sc O + CuS → Sc CuSO
217 -1.28 La SO + CuO → La CuSO
519 -0.54 Sc SO + CuO → Sc CuSO
65 -2.71 LaCuO + 2 La O + La(CuO ) → CuO TABLE IV. Synthesis pathways for various cuprate oxysul-fides based on substitution of sulfur for both (top block) oronly one (middle block) of the apical oxygens in R CuO .Energies in meV/atom and volumes in ˚ A /atom. Since theenergies of formation (∆ E = E products − E reactants ) are pos-itive, none of these pathways appear favorable at ambientconditions. However, high-pressure synthesis will help stabi-lize these pathways, since the majority of volume differentials(∆ V = V products − V reactants ) are negative. We benchmark ourmethod against the standard synthesis pathway for La CuO ,shown on the last line. Surprisingly, ∆ E is +65 meV/atom, soeither DFT systemmatically overestimates enthalpies (whichmeans the actual enthalpies for our hypothetical compoundsare smaller , in our favor), or we must add a bi-directionaluncertainty of ±
70 meV/atom to the computed enthalpies.Additionally, positional entropy of the apical S in the half-substituted R CuSO compounds should also assist in syn-thesis. O LaCuS L a S L a S LaCu LaCu LaCu S O SO SO Cu O CuOCu O L a S L a S C u S La CuSO La CuS O
1: Cu (SO )
2: CuSO
3: La SO
4: La(CuO )
5: LaCuO
6: LaCuO
7: La O
8: La SO
9: LaCuSO10: Cu S
11: LaSO12: CuS13: La CuS O
14: LaCuS
15: La CuS
16: La S O FIG. 22. Gibbs phase diagram of La-Cu-S-O system. Thetwo proposed compositions are marked by red squares. Theyare both found to be unstable: La CuS O lies on the linebetween CuS (12) and La SO (8), and La CuSO lies on thetriangular facet spanned by Cu, La SO (8), and La SO (6).In experiment, we find that the quaternary phase LaCuSO (9)is preferred, likely because copper prefers the 1+ oxidationstate. E n e r g y a b o v e hu ll ( e V / a t o m ) L a S O L a S O L a O L a S C u O C u O C u S L a C u S O L a C u S O L a C u S O FIG. 23. Energies relative to the convex hull for reactants andproducts observed in experiment in the La-Cu-S-O system.Negative energies indicate stable compounds, while unsta-ble compounds have positive energies. Both La CuSO andLa CuS O are highly unstable, lying over 200 meV/atomabove the hull. The vertical axis is broken to display thelarge stability energy of La O . increase in T c were recently provided by STM studieswhich were used to compare the charge transfer gaps ofCa n +1 Cu n O n Cl and Bi Sr Ca n − Cu n O n +4137 , hencethe ideas proposed in Ref. 126 are definitely worth pur-suing. It is not known however, where exactly the max-imum of T c is obtained, as it is clear that a very smallcharge transfer energy would result in an uncorrelatedmaterial with very low T c . So, while this problem remainsopen, we can already say that the work of Ref. 126 hasserved to provide context for later experimental work.Synthesis performed at synchrotron light sources nowallows for in situ tracking of intermediate products andprovides a deep understanding of the chemical reactionpathways by which a compound forms at different tem-peratures . The group of M. Aronson used this tech-nique to study the system described in this section and investigated the potential synthesis of La CuS O and La CuSO . Their work confirmed that the two com-pounds were unstable (at least at high temperatures).From the convex hull (Fig. 22), the target compoundsare thermodynamically unstable with respect to CuS,La SO , Cu, and La SO and the phase with compo-sition LaCuSO seemed to be quite stable at high tem-peratures. In agreement with theory, the end productsLa SO and La SO were observed, in addition to La-CuSO which was the preferred quaternary compositionin almost all the experimentally analyzed reactions. Thisexperimental study thus validates in detail the currentmaterial design workflow and its probabilistic interpreta-tion.Sulfate apical substitution is not easily realized butthe work of Ref. 126 has served to provide context forlater experimental work , which highlight the difficul-ties in forming the desired oxysulfides. The more flexible5 ∆ pd (eV) ∆E (eV) E g (eV) S z ( µ B /Cu) ∆E − band (eV)La CuO CuS O CuS O CuS O CuS O CuO CuO pd = (cid:15) d − (cid:15) p and crystal field splitting ∆E = (cid:15) d ( x − y ) − (cid:15) dz are obtained from the Wannier method applied to GGA(PBE) nonmagnetic calculations. All of the Cu d , O p , and S p orbitals are considered for the Wannier method. The band gap E g and magnetic moment S z (per Cu atom)are caluclated within GGA(PBE) + U antiferromagnetic calculations. The on-site Coulomb interaction U = 8 eV and Hund’scoupling constant J = 0 .
88 eV are used in the GGA+U calculations . The last column of the table describes the results forthe crystal field splitting downfolding to a two band model describing the d x − y and the d z in order to test the theory ofRef. 130. valence of sulfur which can now reduce copper from itsstarting 2+ oxidation state to bring it to its 1+ state,and thus forms LaCuSO instead of our target compounds,which like La CuO , requires a 2+ copper state . Withthis insight we can now go back to the first step of the de-sign loop to search for better routes to reduce the chargetransfer gap while keeping the valence of Cu close to 2+or turn to other routes to vary the charge transfer energyas described in the next section. B. Hg-based Cuprates
With the experience of the previous section we revisitthe question of how to design novel cuprates followingRef. 140. The challenge in the design process is find-ing a new chemical composition that forms in the de-sired cuprate structure - with slightly different parame-ters, thus enabling the elucidation of the mechanism ofsuperconductivity. We focus on the charge transfer en-ergy . We also consider the d z admixture which is thekey variable in the orbital distillation theory . Thistheory is supported by weak coupling calculations thatshow that a reduction in the content of d z (which cor-relates with an increase of the crystal field splitting ∆ E )increases the superconducting critical temperature.Design of novel cuprates is challenging because thephase space is well-explored: most simple point substitu-tions likely have been attempted. On the other hand, anexhaustive search consisting of structure prediction overall possible compositions containing copper would be pro-hibitively costly. Rather than designing novel cupratesfrom scratch, we took an intermediate route: begin withthe family with the highest known transition tempera-tures, the Hg-based cuprates, and modulate its spacerlayers . Setting the Targets, Framing the Questions, HeuristicConsiderations – Taking the Hg-based cuprates as an ex-ample (Fig. 24), we view the cuprates as a stack of func-tional layers, with the composition of each layer chosen to play a specific role. The central copper oxide (CuO )plane supports superconductivity and roughly constrainsthe in-plane lattice constant. The remaining layers musttune the chemical potential of the CuO layer withoutrumpling the plane or introducing disorder, and isolateeach CuO plane to create a 2D system.Our goal is to tune the in-plane charge transfer en-ergy (effective U ), so the relevant layers to focus on arethe BaO layers immediately adjacent to the CuO plane.Due to their spatial proximity, the BaO layers tune thehoppings and interaction strengths of the in-plane Hamil-tonian. We also pay attention to the energy of the d z orbital, as it plays a key role in the orbital distillationtheory .Designing compounds with novel adjacent layers pro-vides a mechanism for controlling superconductivity.However, cycling through all roughly 100 ×
100 elementalsubstitutions for BaO in the periodic table using struc-tural prediction is clearly too naive and computationallyexpensive. To select plausible compositions, we notedthat the BaO layers form a rock salt structure. Usingmaterials databases, we selected all naturally occurringrock salt compounds AX, composed of a cation A and ananion X, starting with 333 in total. We then quickly pre-screened candidates by discarding compositions with (1)large lattice mismatches relative to the in-plane Cu-Cudistance, which we took to be 3.82 ˚A, and (2) anions lesselectronegative than Cu, as these anions would capturedopants intended for the superconducting plane, produc-ing additional Fermi surfaces.
Electronic structure – To evaluate the prospects of su-perconductivity, we follow the previous section, and wefocus on the charge transfer gap. This is summarized inadditional entries in Table V. We find the charge-transferenergy of Hg(CaS) CuO (HCSCO) to be 2.06 eV. It isslightly larger than other Hg-based cuprate, HgBa CuO (HBCO), however is smaller than La CuO and the pro-posed cuprates R CuO S in Table V. Hence we expect T c,max of HCSCO is smaller than that of HBCO, but islarger than those of La CuO and R CuS O (see the6 CuO BaOBaOHgO δ charge dopants structure hamiltonianbalances -2 charge supplies harbors dopants tunes chemical potentialneutral inert protects CuO from disorder tunes in-plane t , t’ , U -2 charge/u.c. accepts roughly sets lattice const. superconducts(same as other CaS layer) FIG. 24. The cuprates are a heterostructure of functional layers, each performing a specific structural and electronic rolein tuning the superconducting hamiltonian. Here, we show the example of HgBa CuO δ , the single layer cuprate with thehighest transition temperature. In our work, we focus on tuning the chemistry of the layers immediately adjacent to the CuO plane, as these will most strongly affect the in-plane Hamiltonian. Adapted from Ref. 140. bottom of Fig. 17).It is also useful to plot the orbitally-resolved bandstructure (Fig. 26). There is indeed a single band crossingthe Fermi level in HCSCO, similar to the other cuprates. Structural Stability – We tested to check that the de-sired structure was stable by first point-substituting theproposed elements into the HgBa CuO (HBCO) struc-ture and checked for local stability via phonon calcula-tions at the Γ, ( π,
0) and ( π, π ) points. If the compositionwas stable in the HBCO structure, then we used USPEXto perform the roughly week-long calculations necessaryto determine whether the HBCO structure is indeed thepreferred energetic minimum.We found three BaO substitutions to be stable in theHBCO structure after phonon screening: CaS, ZrAs, andYbS. USPEX found that only Hg(CaS) CuO (HCSCO)was energetically stable in the layered cuprate structure. Global stability – We construct the convex hull by com-puting the energies of all known compounds in the Hg-Ca-Cu-S-O system, which would produce a 4-dimensionaltetrahedron. In reality, synthesis is performed in an oxy-gen environment parameterized by the chemical potential µ (O ). We find that for all values of µ (O ), HCSCO isunstable, so we pick the value for which the compound isclosest to the convex hull and plot the results in Fig. 25.Prior to examining the global stability of HCSCO,we benchmarked the methodology on HgBa CuO it-self. We find that at fixed oxygen stoichiometry, HBCOlies 74 meV/atom above the hull. Upon modeling thedoped compound using via a 3 × CuO , the metastable structure theoreticallypredicted to be 65 meV/atom above the hull does indeedexist.We found that HCSCO lies 170 meV/atom above thehull (existence probability: 0.09) at fixed oxygen stoi-chiometry, and 240 meV/atom above the hull (existenceprobability: 0.05) for fixed oxygen chemical potential ,which means the compound is likely unstable, but notout of the realm of possibility for successful synthesis.Revising the phase stability calculation - using a moreaccurate correlated method, as suggested in section V, we Cu CaOS OHg HgS Cu S CuS CaSCuS Hg(CaS) CuO FIG. 25. Convex hull for the Hg-Ca-Cu-S chemical system atan oxygen chemical potential of µ (O ) = − .
48 eV, chosenbecause HCSCO is least unstable at this value. The phasediagram forms a tetrahedron with S O, Hg, Cu and CaO atthe vertices (elemental sulfur and calcium are not stable underthis oxygen environment). HCSCO lies in the interior of thetetrahedron, on the triangular face formed by Hg, Cu andCaS. consider the determinant reaction. This reaction is foundby GGA to be HgO + 2 CaS + CuO → Hg(CaS) CuO .Using the GGA(PBE) + Gutzwiller method to calculatethe energies of the reactants, we found that HCSCO isstable with 151 meV/atom below the hull (at fixed oxy-gen stoichiometry). Conclusions – This work in modulating the spacerlayers in the Hg-based cuprates highlights the chal-lenges in optimizing properties in highly-explored ma-terials classes. In the prediction of structures for novelphases, chemical intuition is still crucial for filtering pos-sible candidates and focusing on the most promising com-positions. Structure prediction is the bottleneck step in7 E n e r gy ( e V ) Г ГX M
FIG. 26. Computed band structure of the proposed com-pound Hg(CaS) CuO using GGA(PBE) nonmagnetic calcu-lation. Cu 3 d z and 3 d x − y orbital characters are weightedby red and blue filled-circles, respectively. Similar to thecuprates, a single Cu 3 d x − y band disperses across the Fermilevel. the workflow, consuming the largest fraction of the com-putational resources required. Additionally, known ma-terials can be metastable, and new work must explorewhat factors select for which metastable compounds formin experiment. It would therefore be interesting to pur-sue the MBE route to approach the synthesis of this typeof compound. Further searches which results in similarstructures while avoiding the use of toxic chemicals suchas Hg are also very important.Theoretical calculations reproduced the early observa-tion of Raychaudhury et al. , that one can correlate T c with the oxygen-oxygen overlap and the next-nearest-neighbor hopping parameter of the Hubbard model.Other correlations were pointed out - in particular thatdecreasing the charge transfer gap in these materialsalso increases T c . Finally, the orbital distillation pro-posal argues that a large admixture of the apical or-bitals, in particular the Cu- d z orbital, into the d x − y band, suppresses T c . It has motivated spectroscopicefforts to determine the degree of orbital distillation usingARPES and STM measurements to map the charge-transfer gaps .Finding materials where the admixture of Cu- d z or-bital increases with the smaller charge-transfer gap wouldgreatly advance our understanding of the mechanism ofhigh temperature superconductivity in the cuprates, andas discussed in section VI A, it remains an outstandingchallenge in material design.Exploring the charge-transfer energy (Refs. 126, 145)dependence of T c is worth pursuing further . There areseveral outstanding problems. There should be an opti-mal charge transfer energy to enhance T c , as it is clearthat a very small charge transfer energy would result inan uncorrelated metal with very low T c . It is not known t e m p e r a t u r e pressure, dopingordered phase(CD/CDW) superconductivity FIG. 27. Schematic phase diagram of the potassium-dopedBaBiO system, which is representative of correlated mate-rials as a whole. The phase diagram of correlated systemsgenerally contain an ordered “parent” phase. This parentphase can be suppressed as a function of external parameterslike pressure or doping, and superconductivity can arise nearthe phase transition. however, where exactly the maximum of T c is obtained,as it is clear that a very small charge transfer energywould result in an uncorrelated material with very low T c .Furthermore, taking the results of this section togetherwith the previous section, it seems that the reduction ofthe charge transfer energy makes the formation of thecompound more difficult (i.e. increases its distance fromthe convex hull). Once again, what are the limits whichcan be realized chemically, i.e. how much can the chargetransfer gap be reduced while maintaining a stable com-pound - is a another interesting question. C. Valence Disproportionation in Other HighTemperature Superconductors: CsTlCl Motivation. – In our next example, we seek to design aparent compound for a new superconductor by creating acorrelated material with valence disproportionation andstrong electron phonon coupling. The hypothesis guidingthe project, is that upon suppression of the charge/va-lence order by some means (doping, pressure, disorder),pairing will become the dominant instability and a domeof superconductivity will emerge. This guiding princi-ple is ubiquitous in correlated materials: superconduc-tivity generally appears upon suppression of an ordered“parent” phase as a function of experimental tuning pa-rameters (see the prototypical phase diagram shown inFig. 27). For a recent realization of this idea see Ref. 146.There are cases where the superconducting transitiontemperatures turn out to be well above those expectedfrom the Migdal Eliashberg theory with the electronphonon coupling evaluated within LDA. These materialswere dubbed the “other” high temperature superconduc-tors .A concrete realization of these principles is provided8by Ba − x K x BiO , a well-known superconducting systemdiscovered in the 1980s . The parent compoundis BaBiO , a ∼ . . The order parameter inthis parent state is charge disproportionation, with thebismuth ions nominally alternating between the 3+ and5+ valence (in reality the actual charge valence dispro-porationation is much smaller ). Doping the bariumsite with potassium suppresses the structural distortionsalong with the charge disproportionation and gives riseto superconductivity with a transition temperature ofnearly 30 K at the optimal doping.Conventional electronic structure descriptions basedon LDA/GGA fail to describe the insulating characterof the parent compound. More importantly, the DFT es-timates of the electron phonon coupling λ within Migdal-Eliasberg theory give a value of 0.34 in the doped com-pound, too small to account for its superconductivity .Careful examination found that λ is substantially en-hanced relative to its DFT estimate to a value of nearly1 .
0, and that this enhancement is responsible for super-conductivity in doped BaBiO .It was proposed that static correlations similarly en-hance the electron phonon coupling in other materialsproximate to an insulating state, accounting for super-conductivity in systems such as HfNCl, borocarbides andbuckminsterfullerenes. For these materials, the most im-portant type of correlation that must be captured isthe static contribution, and a GW or hybrid DFT cal-culation is therefore necessary to correct the electronicstructure. After these calculations are done, one isleft with a strongly-coupled electron-phonon system with λ ∼
1. This coupling induces a large dynamical self en-ergy, which accounts for the observed anomalous opticalproperties of doped BaBiO at energies below 1 eV.Hence the low energy electron phonon treatment requiresDMFT.
Setting the targets and Framing the questions. Heuris-tic considerations – Having identified a superconductingmechanism, namely (static) correlation enhanced elec-tron phonon coupling it is natural to seek a realizationof this mechanism in a new material. This task was un-dertaken by Z. Yin et al. in Ref. 153 . The parent com-pound would need to exhibit charge disproportionationand the value of the electron-phonon coupling must beunderestimated by LDA/GGA. The heuristic reasoningused to select CsTlCl follows by analogy with BaBiO .It requires an ion which would charge disproportionate.Like bismuth, thallium is known to valence skip, pre-ferring either a 1+ or 3+ valence state. We want alsothe same perovskite structure which requires the A Tl X composition. Balancing the ionic charges in the presenceof the average 2+ charge of the thallium adds the nextconstraint. Taking ionic radii into consideration, one ar-rives at Cs for the A site, and Cl (or F) for the anion X and lead to the proposal of CsTlCl as a valence dis-proportionation compound, which could be turned (byapplication of pressure and doping) into another “other” FIG. 28. Observed perovskite structure of cubic phase ofCsTlCl . The thallium ion disproportionates in to Tl andTl nominally, and the structure displays an associated al-ternating expansion and contraction of the TlCl octahedralcages. Adapted from Ref. 152. CsTl Cl CsCl T l C l T l C l T l C l Cs TlCl C s T l C l FIG. 29. Gibbs phase diagram of the Cs-Tl-Cl chemi-cal system. The target compound CsTlCl is found to lie3 meV/atom above the convex hull. Although energies abovethe hull strictly imply instability, its small value is indis-tinguishable from zero given the systematic uncertaintiesin current DFT total energies (up to 50 meV/atom), andmany compounds are known to be metastable, lying up to100 meV/atom above the hull. The 3 meV/atom result shouldbe interpreted as a green light to proceed with further inves-tigation as the structure is not obviously unstable. high temperature superconductor . Electronic structure – Identification of the magnitudeand nature of the correlations are an important first stepfor a successful electronic structure prediction, which willlink the structure of the material to the desired property.In order to compute the electronic structure, it was as-sumed that the structure of CsTlCl would be the desiredperovskite crystal structure, shown in Fig. 28.9 FIG. 30. CsTlCl single crystals resulting from successful ma-terial design (left). This material exhibits enhanced electron-phono coupling and is close to metalization with pressure(right), as predicted theoretically (from Ref. 152) To establish that this material is another po-tential “other high temperature superconductor” theelectron-phonon coupling was evaluated using both LDAand methods that capture static correlations, namelyscreened hybrid density functional theory using theHSE06 functional and GW . The calculated electron-phonon couplings at 2.2 GPa with 0.35 hole doping/f.u.(using the virtual crystal approximation) are 0.91 and2.32 for LDA and HSE06 functionals, respectively .The GW method was used as a benchmark to deter-mine the value of the HSE06 screening parameter thatbest reproduced the GW band structure. It was alsoveritfied using this methodology that takes into account static correlations that charge disproportionation occursat half filling .The phase diagram of the material was then deter-mined as a function of pressure and doping. It did re-semble the phase diagram of Fig. 27. Furthermore, theelectron-phonon coupling strength λ attains a value of2.32 at 2.2 GPa, and the predicted superconducting tran-sition temperature is 21 K. These results encouraged thepursuit of experimental synthesis. Structural prediction – Structural stability was testedby computing the most important phonon modes withinDFT and verifying that none were imaginary. It wouldbe interesting to investigate whether this structure wouldbe correctly predicted by the methods described in sec-tion VI E, and what additional polymorphs are possiblein this new class of materials.
Global stability – In the original work of Ref. 153,only a few reaction pathways against known binarieswere checked and the authors reported finding exother-mic reactions. We now check the stability of CsTlCl for decomposition against all known elements, binariesand ternaries in the Cs-Tl-Cl chemical system. Us-ing the computed total energies stored in the MaterialsProject database and the phase diagram function-ality provided with pymatgen , we construct the con-vex hull plotted in Fig. 29. We find that CsTlCl lies3 meV/atom above the convex hull, indicating that thematerial is very likely to form. Lessons learned – Orange crystals of CsTlCl were suc-cessfully grown, which indeed adopt a perovskite struc-ture exhibiting strong breathing distortions indicative FIG. 31. Crystal structures of representative Fe based su-perconductors. They have common FeAs layers but differentspacers. of charge disproportionation at the thallium sites as pre-dicted by theory. The optical gap determined by ab-sortion measurments was of the order of 2.1 eV, in closeagreement with the GW theoretical predictions . Thecrystals are shown in Fig. 30Pressure studies showed that as predicted, the com-pound was very close to metallization. So clearly inthis case, theoretical considerations predicted an interest-ing material which was not listed in the ICSD database,and established that the electron phonon coupling is en-hanced over its LDA value due to static correlations.These halide perovskiets are therefore a new arena toexplore strong electron phonon coupling.This compound turned out to be challenging to dope:doping via Cs vacancies, replacement of Cl by O, S, orN, and substitution of Tl by Hg all did not succeed inpushing the material into a metallic state . This de-serves more careful investigation to see if it is the re-sult of phonon induced self-localization or disorder. Thedopablity of a material, and identification of the mostprobably dopants was intensely studied in weakly corre-lated semiconducting materials and is an outstandingopen problem for correlated insulators and semiconduc-tors worth pursuing.Regardless of whether superconductivity will eventu-ally be discovered in these thallium halides, we found thatwhen motivated by a guiding principle, materials designof a correlated material is possible and yields not onlynew materials but valuable insights into the challenges(some unforseen) required to be overcome for successfulend-to-end materials design. Furthermore, another poly-morph phase was also found in the experiment, but sincean exhaustive structural search had not been performed,it is not possible to ascertain if this phase is among thelow lying structures.
D. Fe 112 compounds
As a final example of the use of material design to elu-cidate the mechanism in correlated superconducting ma-terials, we turn our attention into iron-based supercon-ductors. The discovery of the high temperature super-0 (b) (c)
Γ XY M (ARPES) Ca La FeAs (DMFT) Ca La FeAs (DFT) BaFeAs Γ XMX
Γ X M Γ Z R A Z012-1-2 E n e r gy ( e V ) (a) FIG. 32. Electronic structures of Fe 112 compounds. (a) Band structure (top) and Fermi surface (bottom) of the hypotheticalBaFeAs material, computed by DFT. (b) Angle-resolved photoemission spectroscopy (ARPES) k − E map (top) and thetwo-dimensional (2D) contour of ARPES Fermi surface (bottom) of Ca . La . FeAs . (c) Spectral function A ( k, ω ) (top) andFermi surface (bottom) of Ca . La . FeAs computed by DMFT. The in-plane p orbital characters ( p x and p y ) of the metallicAs spacer are highlighted in red in the band dispersions of (a) and (c). The 2D Fermi surfaces in (a), (b), and (c) are obtainedat k z ∼ π/c . (b) and (c) are adapted from Ref. 158. conducting iron-based materials by Hosono et al. trig-gered a heated debate about the superconducting mech-anism in these materials. The original high T c super-conductor was in the 1111 structure, followed by the 111and the 122 structure (see Fig. 31). In all these cases thespacer layer is inert (LaO for the 1111, Sr++ or Ca++for the 122 case, Li+ in the 111 case) in the sense thattheir electrons are far from the Fermi level.Following the discovery of these materials, multipleideas where put forward in order to understand the ori-gin of their superconductivity. One set of ideas stressesthe importance of the magnetism of the iron ions, as inspin fluctuation theories .A different school of thought posits that what is impor-tant is the high electronic polarizability of the pnictidesand chalcogenides. In the latter case one could envisionthat presence of additional metallic spacer layers wouldmodify the polarizability and will strongly modify thesuperconducting transition temperature as in the pola-ronic mechanism of Ref. 164. This may also be the case ifthe superconductivity is mediated by orbital fluctuations . Designing an iron pnictide superconductor, with ametallic spacer layer and studying how it affects the crit-ical temperature would help elucidate the mechanism ofsuperconductivity in this important class of compounds.This task was undertaken by Shim et al. . Electronic structure – A new family of iron basedsuperconductors, the 112 family was proposed theoret-ically based on chemical analogies and density func-tional studies . Analogies with known 1111 compoundsand the spin fluctuation mechanism lead the authors ofRef. 166 to suggest that BaFeSb and BaFeAs would behigh temperature superconductors, but unlike the knownstructures at the time the 112 structure would have ac- tive spacer layers, where electrons living in the spacerlayers contribute their own Fermi surface as shown inFig. 32(a).Attempts to synthesize iron pnictide materials in the112 structure proposed in Ref. 166 were not originallysuccessful, but new Mn-based materials in this struc-ture were found and both theories and ex-periments show that the spacer layers in theMn-based materials possess Dirac cones. New Ag-basedmaterials in this structure were also reported to possessDirac cones . Structure prediction – The original calculations ofRef. 166 relaxed lattice parameters but did not exam-ine the local stability of the CaAs spacer layer, andhad not considered the stability against phase separa-tion. Motivated by the earlier theoretical work and ex-perimental developments these issues were reexam-ined in Ref. 178, which investigated the crystal struc-ture of CaFeAs by using USPEX in conjunction withVASP as DFT engine. A GGA(PBE) functional and adense Monkhorst-Pack sampling grid with a resolution of2 π × . A − for the k -space integrations. In the struc-tural search in USPEX, two formula-units per unit-cellwere considered. The evolutionary structure predictionpreformed by USPEX predicts that CaFeAs has mono-clinic P /m structure with the distortion of the CaAslayer into zigzag chains. Global stability – A study of the thermodynamic phasestability of CaFeAs was performed in Ref. 178. It is re-produced in Fig. 33. Taking into account stripe mag-netic order in the predicted monoclinic P /m struc-ture of CaFeAs because the additional symmetry break-ing due to magnetic order (for example, antiferromag-1 As FeCaCa As CaAsCa As CaAs FeAs FeAsCaFe As CaFeAs FIG. 33. Ternary phase diagram for CaFeAs . CaFeAs is put above the convex hull and the energy above hull is13 meV/atom. Adapted from Ref. 178. netic order) is not considered in the structural predic-tion performed by USPEX. Due to the stripe magneticorder, the monoclinic P /m structure is further re-laxed into a triclinic P ¯1 structure and the total energyis lower by 19.50 meV/atom (the magnetic moment is1.95 µ B /Fe). CaFeAs lies 13 meV/atom above the con-vex hull, which corresponds to the existence probabilityof 0.37 as discussed in section V. Doping it with rare-earth ions could help its stability with possible energeticgain or an entropic gain regarding no ordering of therare-earth atoms . Similar hypothetical compoundsSrFeAs and BaFeAs were found to be 24 (existenceprobability: 0.33) and 17 meV/atom (existence probabil-ity: 0.36) above the convex hull, respectively , whichare slightly higher in energy compared to CaFeAs .Revising the phase stability calculation - using a moreaccurate correlated method, as suggested in section V,we consider the determinant reaction. This reaction isfound by GGA to be CaAs + FeAs → CaFeAs . The on-site Coulomb interaction U = 5 eV and Hund’s couplingconstant J = 0.8 eV were used in the Gutzwiller calcu-lations . We found that CaFeAs lies 7 meV/atomabove the hull. Electronic structure – The existence of an extra Fermisurface with CaAs character in Ca − x La x FeAs was con-firmed by angle-resolved photoemission and the photoe-mission spectra is in good agreement with DFT+DMFTcalculation (see Figs. 32(b) and (c)). Since the CaAslayer possesses the conducting zig-zag As chain, it in-duces the large electronic anisotropy . Conclusion – Recently, iron-based superconductors inthe 112 structure were synthesized with rare-earth dop-ing, Ca − x La x FeAs and (Ca,Pr)FeAs . Surpris-ingly, these materials form in a structure where the Asin the CaAs layers are distorted in zigzag chains (see
Anion height from Fe layer (Å) Sc O Fe P LaFeAsO F BaFe As P SrFe As (HP)CeFeAsO F NdFeAsO
SmFeAsO F NdFeAsO
TbFeAsO Ba K Fe As Sr V O Fe As FeSe (4 GPa)
BaFe As (HP)NaFeAsLiFeAsFeTe Se FeTe S FeSe (0 GPa) Ca La FeAs Ca Pr FeAs T c on s e t ( K ) Ca La Fe(As Sb ) FIG. 34. T c vs anion height in iron pnictide and chalcogenidesuperconductors. Adapted from Ref. 178. Fig. 31). The space group is monoclinic either P or P /m rather than the originally assumed tetragonalstructure. Second harmonic generation experiments con-firmed the space group P for La-doped compounds ,but similar data are absent for other rare-earth dopingcompounds to the best of our knowledge.Armed with this information, the critical temperatureof the newly discovered 112 compounds was replotted onthe graph of Mizuguchi et al. which displays the T c as a function of the pnictogen height (see Fig. 34). Thepoints do not deviate much from their universal plot, in-dicating that the fate of the superconductivity residesprimarily in the FeAs layers, and is not affected by thenature of the spacer layers. This observation helps ruleout a charge fluctuation mechanism in favor of spin me-diated superconductivity in the iron pnictides.More recently theoretical studies focused on the spacerlayers found that the As p x and p y orbitals are responsiblefor the Dirac cones, and that spin-orbit coupling couldnot only open a gap, but also induce topological phaseson these layers. This suggests that the 112 compoundsare prime candidates for proximity induced topologicalsuperconductivity . More generally the 112 struc-ture provides inspiration for combining iron pnicitide lay-ers with layers having non-trivial topological band struc-ture. This area of research is as yet unexplored, and callsfor a new iteration of the material design loop outlinedin Fig. 9.2 FIG. 35. Crystal Structure of the 2-dimensional parent mate-rial, BaCoS . The green, blue, and yellow spheres correspondto Ba, Co, and S atoms, respectively.FIG. 36. BaCoS phase Diagram from Ref. 183, (c) (1999)The Physical Society of Japan. E. BaCoSO
Motivation – Finally we return to the organizing prin-ciple, shown in Fig. 27, that exotic phases are oftenfound upon suppression of a parent ordered phase oflayered quasi-two-dimensional compounds. The materialBaCoS (Fig. 35) is a layered antiferromagnetic Mott in-sulator. The application of pressure does suppress themagnetic order, but rather than finding a new phase atthe critical point, the material simply becomes metallicdown to the lowest temperatures measured as shown inFig. 36.Clearly, the expectation of novel phases near the sup-pression of order is not a sure bet! Other factors must beinvolved. Still this material is quite interesting, as it dis- FIG. 37. Spectral function for BaCoS at T=180 K usingLDA+DMFT code . plays a transition or crossover from a paramagnetic insu-lating phase to a metallic phase at high temperatures andantiferromagnetically ordered state at low temperatures,features reminiscent of V O and N iSe x S − x , hence acomparative study of these systems can illuminate thefactors that govern the correlation induced metal to in-sulator transition. Targets and questions – The question we address is howthe physical properties of BaCoS change if we substitutethe large sulfur ion by the smaller oxygen ion in BaCoS to form BaCoSO, in particular how does the gap changeas a result. At first sight one would expect that the largersize of the sulfur ion would result in a larger gap, but theanswer to this question is not obvious since structuralchanges can take place. We approached this problemas a case study of our methodology before this materialwas reported in the experimental literature . Some ofthe results were very surprising and are reported in thissection. Electronic Structure – Previous LDA+U and DMFTstudies of BaCoS were reported in Refs. 185 and 186.We returned to this problem using LDA+DMFT usingthe implementation of K. Haule and display the re-sults in Fig. 37. Using the experimental lattice parame-ters, BaCoS is a small gap Mott insulator very close tothe metal insulator transition. Structure prediction – We continue the material de-sign workflow by turning to structural properties of theBa-Co-S-O system by first using USPEX, a structure pre-diction package based on a genetic-search algorithm, tosample the local minima in the energy landscape. Wechoose a 1:1:1:1 ratio for the elements and allowed for twoformula units in a unit cell. We use VASP as our DFT en-gine. We use spin-polarized PBEsol and do not include U corrections. In order not to miss crucial seed structures,which ultimately led to the experimental structure, wefound that the randomly generated initial population ofstructures must be sufficiently large. An initial popula-3 U J (eV)10050050100150200 EE U ( d E / d U ) ( m e V / a t o m ) FIG. 38. The U -dependence of the total energies for 58of the lowest energy structures produced by USPEX dur-ing structure prediction for BaCOSO. The energies are pro-duced by GGA+U allowing for spin polarization. The exper-imentally observed structure (lowest bold purple line) doesnot initially begin with the lowest energy at U = 0, butrapidly becomes favored energetically with increasing correla-tion strength. The group of structures with the lowest ener-gies at U = 0 (yellow lines) rapidly rise nearly linearly with U .Lines colors are graded according to initial slopes at U = 0.The intercept E = − .
92 eV and slope dE /dU = 0 .
19 ofthe energy curve for the experimental structure have beensubtracted for all curves to improve clarity. tion of size 300 was sufficient with a single generationsize of 60. In total ∼
700 metastable structures wereproduced in 8 generations.With GGA the observed structure does not have thelowest energy: it is 271 meV above the lowest-energystructure, and there are 23 structures with lower energy.With spin-polarized GGA the situation is improved, butthe structure is still not among the 10 lowest energy re-sults. A naive conclusion would be that it is not theground state. One would like to distinguish the observedstructure within the low-lying set, for example within0.5 eV/(unit cell) of the final lowest energy structure.However, there are 58 such reported structures (out ofthe final set of 152 best-fit structure), even after remov-ing similar structures.We notice however that adding correlations to theCobalt atom (in the form of a non-zero U ) resolves thisproblem. The energies of the 58 lower-lying structuresare examined as a function of U ( J = 0, which we plot inFig. 38. Crystal structures are kept fixed and only elec-tronic convergence is performed as U is tuned. The linesare colored by their slope at U = 0 which helps visualidentification of the different types of curves.We find that the inclusion of even a very small U ∼ FIG. 39. Experimental structure (left), and representative ofthe group of lowest energy structures found by GGA (right)for the BaCoSO composition. The experimental structurecontains corner-shared CoS O tetrahedra slightly separatedby Ba spacer layers. The GGA structure contains edge-sharing CoS tetrahedra in a fluorite layer, which are sep-arated by BaO rock salt layers. . O tetrahedra forming a cor-rugated 2 D layer, and has been confirmed to be the cor-rect “ground state” structure in experiment. The energygap between the next-best structure and the ground statewidens significantly as U increases. The group of struc-tures the lowest energies initially (set of diagonal yellowlines in Fig. 38) is found to be penalized most rapidlyby correlations, rising nearly linearly. They correspondto slight variations of structures containing CoS fluoritelayers separated by BaO rock salt layers, and a represen-tative is shown in Fig. 39 (right). The remaining struc-tures exhibit intermediate behavior as a function of U .This grouping is most clearly visualized in a plot of theslope vs. the total energy, as shown in Fig. 40. The set ofstructures favored by GGA form a clear outlying clusterin the upper left (yellow), the experimental structure hasthe smallest slope and a competitive total energy (boldsquare, lower left), and the remaining structures are scat-tered in the upper right quadrant.In order to understand this behavior, we examine thedensity matrix n αβσ of cobalt extracted from the DFTcomputations. The indices α and β run over the 3 d or-bitals. The Coulomb U term in LDA+U has the form U − J (cid:80) σ tr { n σ − n σ } , where the trace is over the orbitalindices. The closer the density matrix n σ is to idempo-tency, the less correlations will penalize the state ener-getically.From a physical viewpoint, LDA+U penalizes struc-tures in which the electrons are itinerant. The aver-age occupancy of (spin-resolved) orbitals in itinerant sys-tems cannot be integer: strong charge fluctuations gener-ated by hopping terms prefer occupancies near 0.5. TheLDA+U term penalizes large values of n σ − n σ , whichis maximal when the eigenvalues of n σ ∼ .
5, and thus4 d E / d U FIG. 40. Scatter plot of initial slope of the GGA+U energycurve vs. the total energy at U = 0 for 58 candidate struc-tures produced during the search for the crystal structure ofBaCoSO. The experimental structure (lower left bold purplesquare) does not have the lowest energy, but has the smallestslope, so it is least penalized by increasing U . The group ofstructures with the lowest initial energies (upper left clusterof yellow dots) have large slopes, and are heavily penalizedby correlations. Dotted lines centered on the experimentalstructure are guides to eye. Dots are colored by initial slope,matching the convention in Fig. 40. prefers systems with nearly integer occupancies. We haveplotted the initial slope of the E ( U ) curve against thisLDA+U term in Fig. 41. As expected, the slope andLDA+U term track each other nearly exactly. The ex-perimental structure has orbitals with occupancies near-est to integral, and is least penalized by correlations. Thestructures favored by GGA have the most itinerant or-bitals, likely due to the strong covalent network of CoS tetrahedra in its fluorite block, and are strongly disfa-vored by U .Our proposed strategy for structural prediction is asfollows: first perform USPEX runs with spin-polarizedDFT to generate the list of structures occupying localminima in the energy landscape. Then, apply an accu-rate and likely more computationally expensive method(such as LDA+U or GW) to the resulting structures toreorder the total energies to determine the true groundstate structure. This is more economical than runningUSPEX with hundreds of calls to LDA+U, producingwhat we estimate to be a factor of 5 speed up or more.Finally, we verify that our conclusion about the im-portance of correlation is independent of the methodused. Indeed, as can be seen in Tables VI and VII,the ordering between BaCoSO-expt (the observed struc-ture) and BaCoSO-GGA (the lowest energy structurefound by GGA) is the same in Wien2K as it is in VASPGGA. Calculating the correlated energy in Gutzwiller combined with full-potential linearized augmented planewave Wien2k (W2K) with U = 10 eV flips the order,and makes BaCoSO-expt favorable. Global stability – The convex hull for the Ba-Co-S- n n }1.41.51.61.71.8 d E / d U FIG. 41. Slope of energy vs. U curve plotted against theLDA+U correlation term. Experimental structure is the boldpurple square in the lower left. Dotted line is the diagonal.Dots are colored by initial slope, matching the convention inFig. 40.
O system is shown in Fig. 42. We find BaCoSO lies102 meV/atom above the hull (in the corrected Materials-Project scheme). Observing the importance of correla-tions in determining the correct sign or energy differencesin this system, we examine the effect of correlations onthe convex hull of energies.The relevant reaction for determining the phase stabil-ity of BaCoSO is4 Co + 3 BaS + BaSO → , (31)where left-hand side compounds are on the corner of thegreen triangle in Fig. 42.The W2K and Gutzwiller total energies for materialsin Eq. (31) are listed in Table VII. Note that we used thesame on-site U and Hund’s coupling J in Co d orbitalsfor BaCoSO and Co materials to avoid the GGA/GGA + U correction described in Appendix A. As can be seen,GGA+U as well as ferromagnetic-Gutzwiller calculationsconsistently stabilize the observed material BaCoSO-exptmore than BaCoSO-GGA. Electronic Structure – Knowing the structure we cannow iterate the material design loop and return to thestudy the electronic structure. For this purpose we usethe LDA+DMFT method described in section III. Weshow the calculations of the spectral functions of BaCoSOin Fig. 43 using LDA+DMFT code and compare itto Fig. 37. The combined effects of the substitution andthe structural change result in a Mott insulator with anincreased gap relative to BaCoS . Conclusion – As mentioned earlier, BaCoSO has been5
Ba BaO BaO BaS BaS BaS SCoS Co S Co S CoO CoO Co O CoOS O BaCoSO
1: Ba CoO
2: BaCoO
3: BaCoO BaSO
4: CoSO
5: Co (SO )
6: SO
7: SO FIG. 42. Gibbs phase diagram of Ba-Co-S-O chemical system.The newly synthesized compound, BaCoSO is found to bemarginally unstable, with an energy 102 meV/atom abovethe hull. The compound lies on the triangular facet spannedby Co, BaS and BaSO .FIG. 43. Spectral function for BaCoSO at room temperatureusing LDA+DMFT code . very recently synthesized. It would be very useful to mea-sure its electronic structure by photoemission and com-pare with the theoretical predictions. Ref. 188 theorizedthat it may exhibit high- T c superconductivity. Clearlythe material design loop can continue to the next itera-tion.Application of our proposed workflow raised severalnovel questions. How does U affect the energy land-scape? In particular, does U simply shift the local min-ima relative to one another, or does it create and destroy GGA GGA+UEnergy bare bare corrected
U JE (BaS) -5.106 -5.106 -5.438 0 0 E (Co) -6.799 -7.110 -7.110 0 0 E (BaSO ) -6.554 -6.554 -7.133 0 0 E (BaCoSO-expt) -5.853 -5.580 -6.390 3.32 0 E (BaCoSO-GGA) -5.909 -5.212 -6.021 3.32 0∆ H f (BaCoSO-expt) 0.219 0.570 0.102 - -∆ H f (BaCoSO-GGA) 0.163 0.938 0.471 - -TABLE VI. Comparison of total energies E and energies offormation H f (both in eV/atom) for compounds in the Ba-Co-S-O system via VASP pseudopotential code. The experi-mental and GGA-determined structures, shown in Fig. 39, arelabeled BaCoSO-expt and BaCoSO-GGA. The bare values ofthe outputs of non-spin polarized GGA calculations are pre-sented in the first column. For the uncorrelated calculations,the second column are the bare values of the outputs of spin-polarized GGA(+U) calculations. These energies were mod-ified according to the method used in the Materials Project,described in Sec. IV, to produce the corrected energies shownin the next column. The values of U and J (in eV) usedare listed in the next two columns. In the bottom two rows,we display the formation energies, which is the result of thereaction 3BaS + 4Co + BaSO → minima? Additionally, when is U necessary for correctreordering of the candidate energies? We expect that U is only necessary for compounds containing atoms withpartially-filled d or f shells or magnetic materials, andthis hypothesis deserves to be investigated. Case Studies - wrap-up
The examples in sections VI A–VI E showcase the ques-tion of “how” and “how well” theory can help guide dis-covery of new strongly correlated materials. We foundthat given the current limitations of existing methods,it is useful to quantify the uncertainties by thinking instatistical terms about the space of materials. This ledto Fig. 14, which described the probability for the for-mation of a compound given its calculated energy withvarious methods. This should assist materials scientistsin the search for new compounds, as it provides a lowerbound for the probability of finding something new in ayet unexplored region of material space. As the methodsfor evaluation of total energies improve, this probabil-ity distribution will approach a step function centeredat zero. While the evaluation of this function with theLDA+G method is not feasible at this point, we have re-calculated the energetics of the most-relevant reactionsof the full LDA convex hull study, and summarized thatin Table VIII.Section VI A focused on tuning the charge transfer en-ergy via S substitution to increase the superconductingcritical temperature in La CuS O cuprates. The prob-ability for this material to exist was estimated to be very6 Energy W2K Gutzwiller
U J
PM PM FM E (BaS) -116,175.236 -116,175.236 - 0 0 E (Co) -37,917.854 -37,914.456 -37,914.634 10.0 1.0 E (BaSO ) -40,090.518 -40,090.518 - 0 0 E (BaCoSO-expt) -68,079.015 -68,078.217 -68,078.525 10.0 1.0 E (BaCoSO-GGA) -68,079.077 -68,078.158 -68,078.169 10.0 1.0∆ H f (BaCoSO-expt) 0.106 0.100 -0.209 - -∆ H f (BaCoSO-GGA) 0.044 0.159 0.147 - -TABLE VII. Comparison of total energies E and energies of formation H f (both in eV/atom) for compounds in the Ba-Co-S-Osystem via Wien2k and Gutzwiller calculations. The experimental and GGA-determined structures, shown in Fig. 39, arelabeled BaCoSO-expt and BaCoSO-GGA. In the Gutzwiller calculations, the reference energies were paramagnetic (non-spinpolarized) Wien2K calculations with U = J = 0. Since Wien2K is an all-electron calculation, the energies are extremely largeand not directly comparable to those produced using VASP in Table VI. Shown in the next two columns are the Gutzwillertotal energies for paramagnetic (PM) and ferromagnetic (FM) orderings, respectively, with the Slater-Condon values of U and J shown in the final two columns (in eV). In the bottom two rows, we display the formation energies, which is the resultof the reaction 3BaS + 4Co + BaSO → small, and indeed in-situ studies showed that it decom-poses into the components predicted by the theory .Pursuing this idea led to Hg(CaS) CuO (Section VI B).It is expected to have a higher superconducting transitiontemperature than La CuO , but less than HgBa CuO .Its probability to exist as a ground state structure withinDFT is 0.09, but when post-processing the results withLDA+G, the material is stable, and therefore has a prob-ability larger than 0.5 to exist.Section VI C focused on finding new materials, wherethe electron-phonon coupling is enhanced by static cor-relations. We estimated the probability of CsTlCl toexist to be close to 0.5 and indeed the new materials(including the Flourine cousin) were successfully synthe-sized. Section VI D focused on the iron pnictide, and theidea to create a new family of compounds, the 112 fam-ily, to introduce more polarizable metallic spacer layersbetween the active iron pnictide layers. The probabilityof CaFeAs to form was 0.37, and indeed members of thisfamily were synthesized.Finally in section VI E we studied the Mott insula-tor BaCoSO. We showed that for this material, consid-erations based on LDA or GGA would have been mis-leading both for its ground state structure and stabilityagainst phase separation, but treating dynamical corre-lations gives greatly improved results. Its probability toexist according to the considerations in section V is 0.17,but when the results are post-processed with LDA+G thematerial becomes stable (see Table VIII) and its proba-bility to exist is bigger than 0.5. Indeed this materialwas reported while this paper was in the process of be-ing completed , and was also successfully synthesizedin the lab of M. Aronson .One should stress again the caveats of section V, thatthe considerations leading to these probability estimatesare crude, as they are based on relatively primitive esti-mates of energies, and ignore-finite temperature contri- butions to the entropy - both phononic and electronic.Metastable structures do form, and those are not in-cluded in ground state considerations. Estimating theprobabilities of metastables states to form is problemwhich only now is receiving attention, see Ref. 113.While a compound which is energetically favored willeventually form, the rate of formation is entirely con-trolled by kinetics of nucleation and growth. Theory ofthese phenomena is needed, and the 112 family illustratesthis point. The formation of CaFeAs requires some Ladoping, even though this hardly influences the value ofthe total energy .One can use of the bounds on probability of formationdiscussed in this article conservatively - as a rejectioncriterion. There is little chance for a material to form ifits probability to exist is less than 0.08, as shown clearlyin the examples of La CuS O and La CuSO in boththeory and experiment , while all materials with prob-ability of the order of 0.5 did form. In intermediate casesthe experimentalist needs to weigh the interest in prop-erties of the proposed material against the probability ofits formation.Theory can of course be refined to assist further thesynthesis effort. Knowing that a target compound is un-stable against decomposition, one could raise the chemi-cal potential of the products in those reactions to stabilizethe desired compounds. One can also check if pressurewould stabilize a desirable composition, and whether thecompound would remain metastable once it is synthe-sized under pressure and the pressure is released .Overall, in the five examples discussed in this article,we see a clear consistency between the probabilistic esti-mates and the outcomes of the materials search process.In the process we gained a deeper understanding of howto control charge transfer energies, and how this can beused to test the mechanism of cuprate superconductiv-ity (sections VI A, VI B). We found new compounds on7 Material Determinant reaction ∆ H Materials proj. (eV/atom) P ( x ) ∆ H Gutzwiller (eV/atom)La CuS O La SO + CuS → La CuS O CuSO SO + 4 Cu + La SO → CuSO CuO HgO + 2 CaS + CuO → Hg(CaS) CuO TlCl + Cs TlCl → CaAs + FeAs → CaFeAs → which we could study enhanced electron-phonon couplingby static correlations and its connection to valence dis-proportionation (section VI C), a new family of iron pnic-tide superconductors which can be studied to understandthe mechanism of superconductivity (section VI D), anda new playground to test LDA+DMFT and probe thehandles that can be used to move materials around theMott transition point (section VI E).The examples illustrated how the different method-ological aspects of the treatment of correlations were usedin practice in the different stages of the material designworkflow: (i) heuristic and intuitive considerations, (ii)structure to property relations, (iii) free energy evalua-tions and electronic structure calculations. Structure to property . It is always useful to knowwhether the correlations are static or dynamic, and towhich extent they are local, as this dictates the method-ology used to go from the presumed structure to theproperties. Section VI C treated compounds where thelong range part of the Coulomb interaction is importantand static correlations play a decisive role. This classof materials, which encompass BaBiO , HfNCl and thenewly designed CsTlCl . While within LDA the electron-phonon coupling constant λ was very small, static corre-lations make λ the order of unity, which results in valencedisproportionation as well as high temperature supercon-ductivity with suitable doping. These were treated usingan LAPW implementation of QPGW on the Matsub-ara axis and hybrid DFT functionals . Here there isa clear path for improving the treatment of materialsand estimate the size of the corrections beyond QPGW.Vertex corrected treatments increase the accuracy at anincreased computational cost .Another important direction is to simplify the GW method so as to accelerate its speed - methods such asthe ones proposed in Refs. 193 and 194 can be useful inthis context. In Refs. 53 and 195, which led to the under-standing of the electronic structure of the Ba − x K x BiO system, the HSE screened hybrid functional was used ex-tensively with a range parameter determined by fittingthe results for the more accurate QPGW. This strategywas essential in the design of the perovskite halides de-scribed in section VI C.Sections VI E and VI D illustrated the use of photoe- mission as a theoretical spectroscopy tool to examine thepotential properties of a Mott insulator and a Hund’smetal material. These materials are characterized bystrong dynamical correlations of quite different nature,charge blocking in VI E and spin blocking in VI D. Treat-ing these classes of materials become increasingly difficultwhen we probe longer ranges, lower energies and lowertemperatures, as unexpected emergent phenomena, suchas superconductivity, appear. These more complex func-tionalities at low temperatures - especially predictions ofsuperconducting transition temperatures and other bro-ken symmetry states - are very challenging, in particularin a real materials setting. LDA + cluster DMFT on aplaquette with a zero temperature exact diagonalizationsolver, was the main tool used to guide the project de-scribed in sections VI A - VI B, as it is difficult to reachthe very low temperature regime with Monte-Carlo meth-ods. Work on Monte-Carlo methods to alleviate the mi-nus sign problem and on alternative accurate finite tem-perature impurity solvers to cover this important region,are under active development in many groups around theworld, and advances in this area will have substantial im-pact in searching for interesting low temperature proper-ties of compounds. Structural and Thermodynamical stability – Evalua-tion of free energies is fundamental to structure predic-tion and thermodynamic stability. In the projects de-scribed in this review article only the zero temperatureelectronic energy was estimated. Structural stability re-quires entropic contributions of both vibrational andelectronic origin. For the latter, LDA+DMFT methodswill be very valuable . Since LDA+DMFT calcu-lations of free energies are very time consuming, moreapproximate methods such as LDA+G have been devel-oped . Extension of these methods to finite tem-peratures is an active area of research .Strong correlations are known to affect structural en-ergy differences. A striking example is its influence onthe phase diagram of elemental Plutonium . Other ex-amples abound in transition metal oxides, where im-provements over LDA predictions were made using hy-brid functionals, LDA+U methods , the random phaseapproximation (RPA) , and Monte-Carlo methods ,and Gutzwiller methods .8In section VI E we showed that correlations are impor-tant for structural predictions in a material near the Motttransition - BaCoSO. They can be incorporated first atthe LDA+U or the LDA+G level. While incorporatingthese into every step of an exhaustive search such as US-PEX is expensive, it is enough to consider them on asmaller subset, by post-processing the results of exten-sive LDA studies, as we demonstrated in section VI E.For thermodynamic phase stability, the need to incor-porate correlations is more widely accepted and GGA+Uis used in all the Materials Projects . However, this isnot enough to achieve the necessary accuracy, and empir-ical corrections, which can be significant in magnitude,need to be added.It is important to achieve accuracy using truly ab-initio methods. There is intensive work on using more elab-orate density functionals , or RPA . LDA+DMFTwith U determined from first principles and GW+DMFTare promising directions. This is important, as the val-ues of U currently used for total energy within LDA+Uare quite different from the ones needed to describe thephotoemission spectra in strongly correlated materials.We advocated recomputing total energies with moreadvanced methods as post-processing - after full LDArelaxations and searches. We found that this post-processing to treat correlation, using LDA+U andLDA+G improves bolth structural prediction and theestimation of the probability of formation of this com-pound. In the future, development of forces for LDA+G,improvement in speed in the forces in LDA+DMFT andautomation of user input for electronic methods, couldenable a more self-consistent treatment of correlations instructural relaxation and searches. Developing heuristics and simplified models – In ad-dition to theory and computational tools, rapid mate-rial design requires intermediate layers of inference, ly-ing roughly in the space currently occupied by rules ofthumb such as Pauling’s rules or the Hume-Rothery rulesfor structurally complex alloy phases. One can hope thatthe development of theory and computational approacheswill systematize these rules as well as produce new ones.It is important to gain some understanding of the land-scape in the space of materials, and this is one of themain objectives of the materials genome initiative .This space, as a matter of principle, is infinite as onecan synthesize an infinite number of solids with molec-ular beam epitaxy . Mapping a given material onto amodel Hamiltonian is akin to providing some set of lo-cal coordinates relevant in this space. The hope in thisexercise is that the parametrization captures importantphysics.The projects described in sections VI A–VI B illus-trated the usefulness of mapping onto a model Hamil-tonian, parametrized by t pd , t pp , U, ε d − ε p and how theycan be used to explore the factors that control T c . Justlike in the case of the iron pnictides (Fig. 34), furtherfundamental understanding of the mechanism is needed.We have argued that designing and synthesizing new ma- terials will help in this process. VII. CONCLUSION AND OUTLOOK
Theory and computation have progressed to the pointthat they are now set to play an important role in mate-rial design and exploration, at least for weakly correlatedelectron materials. For strongly correlated materials, weargued that the material design process is intellectuallyimportant as it serves as a strong test of physical or chem-ical intuition, the quality of our methodology and predic-tive capabilities. We also showed through an admittedlyvery small number of concrete examples, that existingmethodologies already have a potential to accelerate ma-terial discovery.The materials considered in the five exemplars wererelatively small variations with respect to known com-pounds. In this context LDA+DMFT has substantialpredictive power, as one can use a fixed value of the in-teraction parameters, applied to a localized orbital sup-ported on a large energy window. The LDA+DMFTdouble counting correction can then be determined fullyfrom first principles as proposed in Ref. 51. Improvingthe accuracy of the estimates of the remaining parame-ters such as the strength of the screened Hubbard U andHund’s coupling J and the rest of the Slater parametersin DMFT from ab-initio would be very important formaterial design as we could contemplate property pre-diction in completely unexplored regions in the space ofmaterials.Combining multiple electronic structure methods totreat correlations constitutes a major software develop-ment challenge. Advances in this area will facilitate thedevelopment of quantitative measures in the space of ma-terials, and streamline the study of materials. Such toolswould enable to scale investigations, such as the onespresented in this review, from five exemplars to a sta-tistically significant sample. It is also desirable to en-able integration between electronic structure tools andthe growing databases of materials and their properties.More standardization of data and databases will lead toeasier integration between software and databases. Ef-forts in this direction have been initiated . Finally,machine learning methods will provide alternative routesfor exploring structure and property relations.In solid state physics there is a long tradition of fruitfulinteractions between experiment and theory which hasresulted in remarkable advances in understanding con-densed matter systems. Until very recently however, ma-terials discovery was entirely driven by experiment. Thissituation is rapidly changing. The highest temperaturesuperconductor, HS under pressure, was synthesized following a material specific theoretical suggestion .When theory-assisted material design becomes reality,theory will be able to play new additional roles. Theorycan help guide experiment in selecting which composi-tions among the thousands of possibilities would have9the highest probability of forming a new material, andthus be fast-tracked for exploration over less promisingcompositions. In return, experiment can provide tests oftheoretical predictions of structure and properties, whilenew advances in experimental techniques, where the in-termediate products can be monitored in real time withinthe reaction vessel, will shed light into how material syn-thesis actually takes place. This will be a golden agefor experiment - theory interactions, which has potentialto take the field of correlated-electron systems to a newlevel. VIII. ACKNOWLEDGEMENTS
We are grateful to our theory (Zhiping Yin, TuranBirol, C. Weber and Jihoon Shim) and experimental(Martha Greenblatt, Meigan Aronson, Ni Ni, CedomirPetrovic and Emilia Morosan) colleagues for numerouscollaborations which shaped our understanding of thematerial design challenge. We are grateful to Karin Rabe,Kathi Stadler, Piers Coleman and Vladan Stevanovic foruseful comments on the manuscript. This work was su-ported by the US Department of Energy, Office of ScienceBasic Energy Science as a part of the Computational Ma-terials Science Program.
Appendices A. EMPIRICAL CORRECTIONS
We present a few of the empirical correction schemesused in literature below. All three schemes that werepresented corrected GGA with GGA+U for certain ele-ments in specific configurations. There is not completeagreement between the methods on which compounds toapply U to, or what its value should be. Whereas FERE,building on top of Lany’s observation, fitted the Uto prefer the correct stoichiometry between oxides withdifferent stoichiometries, Materials Project fitted the U’son binary oxide formation enthlapies, arriving at differentvalues. Whereas FERE applies U’s universally for the setof chosen elements (although it is fitted once on oxides),Materials Project and OQMD apply GGA+U only forcompounds containing oxygen or flouride. The values ofchosen U (or U eff = U − J in Dudarev’s method ) in allscheme range between 3-4eV in all schemes (other thanone value).In the fitted elemental-phase reference en-ergy (FERE) scheme , experimental formation energies∆ H exp are used to determine the best “energies” of theelemental phases. These fitted elemental-phase referenceenergies E FERE are constructed to minimize the system-atic error in the formation energies of the training set ofcompounds, and can then be used to predict the forma-tion energies of new compounds. The procedure requires the tabulation of the experimental formation energies for252 binary compounds A m B n along with GGA+U cal-culations for each of the compounds to obtain the totalenergies E GGA+U . From this data, the elemental energies E FERE of the 50 elements which span the set of binariesare fitted by solving the linear least-squares problem:∆ H exp (A m B n ) = E GGA+U (A m B n ) − mE FERE ( A ) − nE FERE ( B ) . (32)As emphasized by the authors, the energies E FERE arenot meant to improve the absolute total energies for theelemental phases, but rather constructed to optimize thesystematic cancellation of errors with the GGA+U totalenergies. A final detail is the choice of U in the GGA+Ucalculations. As detailed in Ref. 14, three values are used:3 eV for most transition metals, 5 eV for Cu and Ag,and 0 for the remaining elements. The scheme performsreasonably well: the mean absolute error of the binariesis 260 meV/atom when computed using GGA+U, anddrops to 54 meV/atom in the FERE scheme. These im-provements carry over when the fitted elemental energiesare used to compute the formation energies of a test setof 55 ternaries, with the mean absolute error lying at48 meV/atom.In the framework of Open Quantum MaterialsDatabase (OQMD, Ref. 4), energies of 1,670 reactionswere collected, together with the corresponding GGA orGGA+U results. Unlike the FERE scheme , they ap-plied GGA+U only to a selected set of “correlated” ele-ments when they are in oxides (V, Cr, Mn, Fe, Co, Ni,Cu, Th, U, Np, Pu). Fitting elemental chemical po-tential in experimental formation energies, similarly toabove Eq. (32), they compare two different empirically-based correction schemes. In the ‘fit-partial’ scheme,they fit elemental energies only for those elements wherethe DFT energy of the T = 0 K ground state is notan accurate reference for room-temperature and ambi-ent pressure (STP) formation energies. These elementsinclude room-temperature diatomic gases (H, N, O ,F,Cl), room temperature liquids (Br, Hg), molecular solids(P, S, I), elements with phase transformation between 0K and room temperature (Na, Ti, Sn), and the “corre-lated” elements listed above. Their ‘fit-all’ scheme onthe other hand, allows all elemental chemical potentialsto be fitted. Interestingly they also find that there isconsiderable error within the two different experimentaldata sets in their database. It would be interesting tounderstand the source of this experimental discrepancy.In the Materials Project , an elaborate scheme ofphysically-motivated energy shifts are used to respec-tively correct the formation energy of gases, compoundscontaining electronegative anions and energetic contribu-tions due to correlations at transition metal sites .The corrections in this scheme are classified as follows:1.
Gas Correction - formation-energies for gaseous el-ements are used as fitting variables for experimentaldata, similarly to the other approaches. Thereforethe energies for N , F , Cl and H are tabulated.02. Anion Correction – DFT tends to overbind the O molecule. It requires too much energy to dissociateO molecule in oxidation reactions, which resultsin underestimating the oxidation energy. It givesa constant energy shift compared to experimentaldata . Hence, the DFT total energy for any oxideor sulfide compound is corrected by the followingnegative shift, which accounts for DFT overbindingof anions: ∆ E = E · N, (33)where E is the correction energy, and N is thenumber of O or S in the composition. If bothoxygen and sulfur are present in the compound,both corrections are applied. The correction en-ergy E is tabulated for the ionic state of the an-ion: oxide(O − ), peroxide(O − ), superoxide(O − ),ozonide(O − ), and sulfide(S − ).3. Correction for correlations – GGA+U is applied fora selected set of transition metal compounds con-taining oxygen or fluorine atoms (transition metaloxides or transition metal fluorides) . Thevalue of U was fixed empirically for the followingset of elements: Co, Cr, Fe, Mn, Mo, Ni, V, W. 4.
Correction for GGA+U vs. GGA compatibility –The DFT energy for any oxide or fluoride com-pound is corrected if the run is a GGA+U calcu-lation. The correction is applied to any fluorideor oxide containing one of the following transitionmetals: { V, Cr, Mn, Fe, Co, Ni, Mo, W } . NoteTi and Cu are absent. The corrections are of theform: ∆ E = (cid:88) i E i · N i , (34)where i runs over the transition metals, N i is thenumber of atoms of the transition metal present inthe compound, and E i is the correction energy istabulated.While the corrections are semi-empirical, they per-form well from a practical standpoint, allowing forma-tion energies to be accurately determined to within 25-50 meV/atom as compared to experimental thermo-dynamic benchmarks and successfully reproduce Gibbsphase diagrams . In the paper we used the Ma-terials Project software to construct convex-hulls andextract energies-above-hull . Therefore the correc-tions that we apply are based on their scheme. L. H. Greene, H. Z. Arham, C. R. Hunt, and W. K. Park,J Supercond Nov Magn , 2121 (2012). K. Lejaeghere, G. Bihlmayer, T. Bj¨orkman, P. Blaha,S. Bl¨ugel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark,A. Dal Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst,I. Di Marco, C. Draxl, M. Du(cid:32)lak, O. Eriksson, J. A.Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi,M. Giantomassi, S. Goedecker, X. Gonze, O. Gr˚an¨as,E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann,P. J. Hasnip, N. A. W. Holzwarth, D. Iu¸san, D. B.Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik,E. K¨u¸c¨ukbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck,M. Marsman, N. Marzari, U. Nitzsche, L. Nordstr¨om,T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans,M. I. J. Probert, K. Refson, M. Richter, G.-M. Rignanese,S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma,F. Tavazza, P. Thunstr¨om, A. Tkatchenko, M. Torrent,D. Vanderbilt, M. J. van Setten, V. Van Speybroeck, J. M.Wills, J. R. Yates, G.-X. Zhang, and S. Cottenier, Science , aad3000 (2016). A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards,S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, andK. a. Persson, APL Materials , 011002 (2013). S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W.Doak, M. Aykol, S. R¨uhl, and C. Wolverton, npj Com-putational Materials , 15010 (2015). W. Setyawan, R. M. Gaume, S. Lam, R. S. Feigelson, andS. Curtarolo, ACS Comb. Sci. , 382 (2011). D. Neugebauer, J, Kormann F, Grabowski, B, Hickel T,Raabe, in
Harnessing Mater. Genome Accel. Mater. Dev.via Comput. Exp. Tools , 2013 (2012). C. J. Fennie, Phys. Rev. Lett. , 167203 (2008). R. Gautier, X. Zhang, L. Hu, L. Yu, Y. Lin, T. O. L.Sunde, D. Chon, K. R. Poeppelmeier, and A. Zunger,Nature Chemistry , 308 (2015). D. J. Fredeman, P. H. Tobash, M. A. Torrez, J. D. Thomp-son, E. D. Bauer, F. Ronning, W. W. Tipton, S. P. Rudin,and R. G. Hennig, Phys. Rev. B , 224102 (2011). N. P. Armitage, E. J. Mele, and A. Vish-wanath, Reviews of Modern Physics (2018),https://doi.org/10.1103/RevModPhys.90.015001,arXiv:1705.01111. D. Duan, Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao,H. Yu, B. Liu, W. Tian, and T. Cui, Sci. Rep. , 6968(2014). A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Kseno-fontov, and S. I. Shylin, Nature , 73 (2015). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). V. Stevanovi´c, S. Lany, X. Zhang, and A. Zunger, Phys.Rev. B , 115104 (2012). DMFT at 25: Infinite Dimensions , Modeling and Simula-tion, Vol. 4, Autumn School on Correlated Electrons, Jlich(Germany), 15 Sep 2014 - 19 Sep 2014 (Forschungszen-trum Jlich Zentralbibliothek, Verlag, Jlich, 2014). M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett. ,1418 (1985). S. G. Louie, “First-principles theory of electron excita-tion energies in solids, surfaces, and defects,” in
Topicsin Computational Materials Science , edited by C. Fong(World Scientific, Singapore, 1997) pp. 96,. I. Shavitt and R. Bartlett,
Many-Body Methods in Chem-istry and Physics: MBPT and Coupled-Cluster The-ory , Cambridge Molecular Science (Cambridge UniversityPress, 2009). J. M. Tomczak, M. van Schilfgaarde, and G. Kotliar,Phys. Rev. Lett. , 237010 (2012). M. V. Sadovskii, I. A. Nekrasov, E. Z. Kuchinskii, T. Pr-uschke, and V. I. Anisimov, Phys. Rev. B - Condens.Matter Mater. Phys. , 1 (2005), arXiv:0508585 [cond-mat]. L. Hedin, Phys. Rev. , A796 (1965). B. Holm and U. von Barth, Phys. Rev. B , 2108 (1998). A. Stan, N. E. Dahlen, and R. van Leeuwen, Europhys.Lett. , 298 (2006). A. Kutepov, S. Y. Savrasov, and G. Kotliar, Phys. Rev. B , 041103(R) (2009). S. V. Faleev, M. van Schilfgaarde, and T. Kotani, Phys.Rev. Lett. , 126406 (2004). A. Kutepov, K. Haule, S. Y. Savrasov, and G. Kotliar,Phys. Rev. B , 155129 (2012). A. L. Kutepov, V. S. Oudovenko, and G. Kotliar, Com-put. Phys. Commun. , 407 (2017). G. Kotliar and D. Vollhardt, Physics Today , 53 (2004). A. Georges and G. Kotliar, Phys. Rev. B , 6479 (1992). W. Metzner and D. Vollhardt, Phys. Rev. Lett. , 324(1989). E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov,M. Troyer, and P. Werner, Rev. Mod. Phys. , 349(2011). V. I. Anisimov, A. I. Poteryaev, M. A. Korotin, A. O.Anokhin, and G. Kotliar, J. Phys. Condens. Matter. ,7359 (1997). A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B ,6884 (1998). S. Y. Savrasov and G. Kotliar, Phys. Rev. B , 245101(2004). S. Y. Savrasov, G. Kotliar, and E. Abrahams, Nature , 793 (2001). R. Chitra and G. Kotliar, Phys. Rev. B , 115110 (2001). R. Chitra and G. Kotliar, Phys. Rev. B , 12715 (2000). X. Dai, S. Y. Savrasov, G. Kotliar, A. Migliori, H. Led-better, and E. Abrahams, Science , 953 (2003). I. Paul and G. Kotliar, Eur. Phys. J. B , 189 (2006). J. Lee and K. Haule, Phys. Rev. B , 155144 (2015). B. Amadon, F. Lechermann, A. Georges, F. Jollet, T. O.Wehling, and A. I. Lichtenstein, Phys. Rev. B , 205112(2008). B. Amadon, J. Phys. Condens. Matter. , 075604 (2012). H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev.B , 235103 (2014). O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann,I. Krivenko, L. Messio, and P. Seth, Comput. Phys. Com-mun. , 398 (2015). “Amulet: http://amulet-code.org/,” . “Lisa: http://dmft.rutgers.edu/lisa/,” . “Lmtart: http://mindlab.physics.ucdavis.edu/materialresearch/scientific/index lmtart.htm,” . K. Haule, C. H. Yee, and K. Kim, Phys. Rev. B ,195107 (2010). “Rutgers dft+edmft: http://hauleweb.rutgers.edu/tutorials/,”. O. Gr˚an¨as, I. D. Marco, P. Thunstr¨oma, L. Nordstr¨om,O. Eriksson, T. Bj¨orkman, and J. M. Wills, Comput.Mater. Sci. , 295 (2012). K. Haule, Phys. Rev. Lett. , 196403 (2015). C.-O. Almbladh, U. von Barth, and R. van Leeuwen, Int.J. Mod. Phys. B , 535 (1999). Z. P. Yin, A. Kutepov, and G. Kotliar, Phys. Rev. X ,021011 (2013). D. Jacob, K. Haule, and G. Kotliar, EPL , 57009(2008). S. Biermann, F. Aryasetiawan, and A. Georges, Phys.Rev. Lett. , 086402 (2003). G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko,O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. ,865 (2006). P. Sun and G. Kotliar, Phys. Rev. B , 085120 (2002). P. Sun and G. Kotliar, Phys. Rev. Lett. , 196402 (2004). G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin,A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N.Rubtsov, and K. Held, Rev. Mod. Phys. , 025003(2018). S. Choi, A. Kutepov, K. Haule, M. van Schilfgaarde, andG. Kotliar, npj Quantum Materials , 16001 (2016). J. M. Tomczak, M. Casula, T. Miyake, F. Aryasetiawan,and S. Biermann, EPL , 67001 (2012). V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys.Rev. B , 943 (1991). S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J.Humphreys, and A. P. Sutton, Phys. Rev. B , 1505(1998). A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys.Rev. B , R5467(R) (1995). V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein,J. Phys. Condens. Matter , 767 (1997). A. G. Petukhov, I. I. Mazin, L. Chioncei, and A. I. Licht-enstein, Phys. Rev. B , 153106 (2003). M. T. Czy´zyk and G. A. Sawatzky, Phys. Rev. B , 14211(1994). B. Himmetoglu, A. Floris, S. de Gironcoli, andM. Cococcioni, Int. J. Quantum Chem. , 14 (2014),arXiv:1309.3355. G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. ,1362 (1986). T. Li, P. W¨olfle, and P. J. Hirschfeld, Phys. Rev. B ,6817 (1989). R. Fr´esard and P. W¨olfle, Int. J. Mod. Phys. B , 685(1992). R. Fr´esard and G. Kotliar, Phys. Rev. B , 12909 (1997). F. Lechermann, A. Georges, G. Kotliar, and O. Parcollet,Phys. Rev. B , 155102 (2007). R. Raimondi and C. Castellani, Phys. Rev. B ,11453(R) (1993). N. Lanat`a, Y. Yao, C.-Z. Wang, K.-M. Ho, andG. Kotliar, Phys. Rev. X , 011008 (2015). N. Lanat`a, Y. Yao, X. Deng, V. Dobrosavljevi´c, andG. Kotliar, Phys. Rev. Lett. , 126401 (2017). J. B¨unemann and F. Gebhard, Phys. Rev. B , 193104(2007). X. Deng, X. Dai, and Z. Fang, EPL , 37008 (2008). K. M. Ho, J. Schmalian, and C. Zhuang, Phys. Rev. B , 073101 (2008). S. Y. Savrasov, V. Oudovenko, K. Haule, D. Villani, andG. Kotliar, Phys. Rev. B , 115117 (2005). C. Piefke and F. Lechermann, Phys. Rev. B (2017), https://doi.org/10.1103/PhysRevB.97.125154,arXiv:1708.03191. G. Baym and L. P. Kadanoff, Phys. Rev. , 287 (1961). G. Baym, Phys. Rev. , 1391 (1962). F. Krien, E. G. C. P. van Loon, H. Hafermann, J. Otsuki,M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B , 075155 (2017). G. Biroli, O. Parcollet, and G. Kotliar, Phys. Rev. B ,205108 (2004). D. Karlsson and R. van Leeuwen, Phys. Rev. B , 125124(2016). G. Stefanucci, Y. Pavlyukh, A.-M. Uimonen, and R. vanLeeuwen, Phys. Rev. B , 115134 (2014). A. Kutepov, K. Haule, S. Y. Savrasov, and G. Kotliar,Phys. Rev. B , 045105 (2010). N. E. Zein, S. Y. Savrasov, and G. Kotliar, Phys. Rev.Lett. , 226403 (2006). F. Nilsson, L. Boehnke, P. Werner, and F. Aryasetiawan,Phys. Rev. Materials , 043803 (2017). T. Ayral, T.-H. Lee, and G. Kotliar, Phys. Rev. B ,235139 (2017). S. M. Woodley and R. Catlow, Nat. Mater. , 937 (2008). J. Pannetier, J. Bassas-Alsina, J. Rodriguez-Carvajal,and V. Caignaert, Nature , 343 (1990). S. Kirkpatrick, C. D. G. Jr., and M. P. Vecchi, Science , 671 (1983). S. M. Woodley, P. D. Battle, J. D. Gale, and C. R. A.Catlow, Phys. Chem. Chem. Phys. , 2535 (1999). A. R. Oganov and C. W. Glass, J. Chem. Phys. ,244704 (2006). G. Trimarchi and A. Zunger, Phys. Rev. B , 104113(2007). S. Curtarolo, D. Morgan, K. Persson, J. Rodgers, andG. Ceder, Phys. Rev. Lett. , 135503 (2003). C. C. Fischer, K. J. Tibbetts, D. Morgan, and G. Ceder,Nat. Mater. , 641 (2006). A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. , 12562 (2002). A. Barducci, M. Bonomi, and M. Parrinello, Comput.Mol. Sci. , 826 (2011). D. J. Wales and J. P. K. Doye, J. Chem. Phys. A ,5111 (1997).
S. Goedecker, J. Chem. Phys. , 9911 (2004).
C. J. Pickard and R. J. Needs, J. Phys. Condens. Matter. , 053201 (2011). V. Stevanovi´c, Phys. Rev. Lett. , 075503 (2016).
A. O. Lyakhov, A. R. Oganov, H. T. Stokes, and Q. Zhu,Computer Physics Communications , 1172 (2013).
G. Hautier, S. P. Ong, A. Jain, C. J. Moore, andG. Ceder, Phys. Rev. B , 155208 (2012). The SQL query sent to OQMDv1.1 was:selectd.composition id, ’
A. Watson and T. Markus,
Subvolume C Ternary SteelSyst. Phase Diagrams Phase Transit. Data, Part 2Ternary Syst. from Cr-Mn-N to Ni-Si-Ti Ser. Landolt-B¨ornstein Numer. Data Funct. Relationships Sci. Tech-nol. - New Ser. Part 19C2 (2015).
P. Nash, “Thermodynamic database,” (2013).
E. Jones, T. Oliphant, P. Peterson, et al. , “SciPy: Opensource scientific tools for Python,” (2001–).
A. Narayan, A. Bhutani, S. Rubeck, J. N. Eckstein, D. P.Shoemaker, and L. K. Wagner, Phys. Rev. B , 045105(2016). W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain,W. D. Richards, A. C. Gamst, K. A. Persson, andG. Ceder, Science Advances , e1600225 (2016). P. J. Ryan, J. W. Kim, T. Birol, P. Thompson, J. H. Lee,X. Ke, P. S. Normile, E. Karapetrova, P. Schiffer, S. D.Brown, C. J. Fennie, and D. G. Schlom, Nat. Commun. , 1 (2013). L. Havl´ak, J. F´abry, M. Henriques, and M. Duˇsek, ActaCrystallogr. Sect. C Struct. Chem. , 623 (2015). X. Deng, K. Haule, and G. Kotliar, Phys. Rev. Lett. ,176404 (2013).
Q. Yin, A. Kutepov, K. Haule, G. Kotliar, S. Y. Savrasov,and W. E. Pickett, Phys. Rev. B , 195111 (2011). D. Butler, Nature , 238 (2004).
P. S¨oderlind, G. Kotliar, K. Haule, P. M. Oppeneer, andD. Guillaumont, MRS Bull. , 883 (2010). Y. Yun and P. M. Oppeneer, MRS Bull. , 178 (2011). J. Zaanen, G. A. Sawatzky, and J. W. Allen, Phys. Rev.Lett. , 418 (1985). N. E. Bickers, D. J. Scalapino, and R. T. Scalettar, Int.J. Mod. Phys. B , 687 (1987). G. Kotliar and J. Liu, Phys. Rev. B , 5142 (1988). F. C. Zhang, C. Gros, T. M. Rice, and H. Shiba, Su-percond. Sci. Technol. , 36 (1988), arXiv:0311604 [cond-mat]. M. Capone and G. Kotliar, Phys. Rev. B , 1 (2006),arXiv:0603227 [cond-mat]. C. Weber, C. Yee, K. Haule, and G. Kotliar, EPL ,37001 (2012).
E. Pavarini, T. Saha-Dasgupta, O. Jepsen, and O. K.Andersen, Phys. Rev. Lett. (2001), 10.1103/Phys-RevLett.87.047003, arXiv:0406357 [cond-mat]. Y. Ohta, T. Tohyama, and S. Maekawa, Phys. Rev. B , 2968 (1991). R. Raimondi, J. H. Jefferson, and L. F. Feiner, Phys.Rev. B , 8774 (1996). H. Sakakibara, H. Usui, K. Kuroki, R. Arita, and H. Aoki,Phys. Rev. Lett. , 057003 (2010).
N. Marzari and D. Vanderbilt, Phys. Rev. B , 12847(1997). I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B , 035109 (2001). A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Van-derbilt, and N. Marzari, Comput. Phys. Commun. ,685 (2008).
C. H. Yee and G. Kotliar, Phys. Rev. B - Condens. MatterMater. Phys. (2014), 10.1103/PhysRevB.89.094517,arXiv:1305.0322. V. I. Anisimov, M. A. Korotin, A. S. Mylnikova, A. V.Kozhevnikov, D. M. Korotin, and J. Lorenzana, Phys.Rev. B , 172501 (2004). H. He, C.-H. Yee, D. E. McNally, J. W. Simonson, S. Zell-man, M. Klemm, P. Kamenov, G. Geschwind, A. Ze-bro, S. Ghose, J. Bai, E. Dooryhee, G. Kotliar, andM. C. Aronson, Proc. Natl. Acad. Sci. (2018), 10.1073/p-nas.1800284115.
W. Ruan, C. Hu, J. Zhao, P. Cai, Y. Peng, C. Ye, R. Yu,X. Li, Z. Hao, C. Jin, X. Zhou, Z.-Y. Weng, and Y. Wang,Science Bulletin , 1826 (2016). R. Molinder, T. P. Comyn, N. Hondow, J. E. Parkerc,and V. Duponta, Energy Environ. Sci. , 8958 (2012). D. P. Shoemaker, Y.-J. Hu, D. Y. Chung, G. J. Halder,P. J. Chupas, L. Soderholm, J. F. Mitchell, and M. G.Kanatzidis, Proc. Natl. Acad. Sci. USA , 10922(2014).
C.-H. Yee, T. Birol, and G. Kotliar, EPL , 17002(2015).
S. P. Ong, L. Wang, B. Kang, and G. Ceder, Chemistryof Materials , 1798 (2008). M. De Raychaudhury, E. Pavarini, and O. K. Andersen,Phys. Rev. Lett. , 126402 (2007). C. E. Matt, D. Sutter, A. M. Cook, Y. Sassa, M. Mans-son, O. Tjernberg, L. Das, M. Horio, D. Destraz, C. G.Fatuzzo, K. Hauser, M. Shi, M. Kobayashi, V. Strocov,P. Dudin, M. Hoesch, S. Pyon, T. Takayama, H. Tak-agi, O. J. Lipscombe, S. M. Hayden, T. Kurosawa,N. Momono, M. Oda, T. Neupert, and J. Chang, Na-ture Communications , 1 (2018), arXiv:1707.08491. J. A. Slezak, J. Lee, M. Wang, K. McElroy, K. Fujita,B. M. Andersen, P. J. Hirschfeld, H. Eisaki, S. Uchida,and J. C. Davis, Proc. Natl. Acad. Sci. , 3203 (2008).
C. Weber, K. Haule, and G. Kotliar, Phys. Rev. B ,125107 (2010). L. Li, X. Deng, Z. Wang, Y. Liu, M. Abeykoon, E. Doory-hee, A. Tomic, Y. Huang, J. B. Warren, E. S. Bozin,S. J. L. Billinge, Y. Sun, Y. Zhu, G. Kotliar, and C. Petro-vic, npj Quantum Mater. , 11 (2017). W. E. Pickett, Physica B: Condens. Matter. , 112(2001).
A. Sleight, J. Gillson, and P. Bierstedt, Solid State Com-munications , 27 (1975). R. J. Cava, B. Batlogg, J. J. Krajewski, R. Farrow,L. W. Rupp, A. E. White, K. Short, W. F. Peck, andT. Kometani, Nature , 814 (1988).
V. Meregalli and S. Y. Savrasov, Phys. Rev. B , 14453(1998). R. Nourafkan, F. Marsiglio, and G. Kotliar, Phys. Rev.Lett. , 017001 (2012).
M. Retuerto, T. Emge, J. Hadermann, P. W. Stephens,M. R. Li, Z. P. Yin, M. Croft, A. Ignatov, S. J. Zhang,Z. Yuan, C. Jin, J. W. Simonson, M. C. Aronson, A. Pan,D. N. Basov, G. Kotliar, and M. Greenblatt, Chemistryof Materials , 4071 (2013). Z. P. Yin and G. Kotliar, EPL , 27002 (2013).
S. P. Ong, A. Jain, G. Hautier, B. Kang, and G. Ceder,Electrochemistry Communications , 427 (2010). S. P. Ong, W. D. Richards, A. Jain, G. Hautier,M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A.Persson, and G. Ceder, Computational Materials Science , 314 (2013). M. Retuerto, Z. Yin, T. J. Emge, P. W. Stephens, M.-r.Li, T. Sarkar, M. C. Croft, A. Ignatov, Z. Yuan, S. J.Zhang, C. Jin, R. Paria Sena, J. Hadermann, G. Kotliar,and M. Greenblatt, Inorganic Chemistry , 1066 (2015). A. Zunger, Appl. Phys. Lett. , 57 (2003). S. Jiang, C. Liu, H. Cao, T. Birol, J. M. Allred, W. Tian,L. Liu, K. Cho, M. J. Krogstad, J. Ma, K. M. Taddei,M. A. Tanatar, M. Hoesch, R. Prozorov, S. Rosenkranz,Y. J. Uemura, G. Kotliar, and N. Ni, Phys. Rev. B ,054522 (2016). Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono,J. Am. Chem. Soc. , 3296 (2008).
I. I. Mazin, D. J. Singh, M. D. Johannes, and M. H. Du,Phys. Rev. Lett. , 057003 (2008).
K. Kuroki, S. Onari, R. Arita, H. Usui, Y. Tanaka,H. Kontani, and H. Aoki, Phys. Rev. Lett. , 087004(2008).
S. Graser, G. R. Boyd, C. Cao, H.-P. Cheng, P. J.Hirschfeld, and D. J. Scalapino, Phys. Rev. B ,180514(R) (2008). A. V. Chubukov, D. V. Efremov, and I. Eremin, Phys.Rev. B , 134512 (2008). M. Berciu, I. Elfimov, and G. A. Sawatzky, Phys. Rev. B , 214507 (2009). S. Onari and H. Kontani, Phys. Rev. Lett. , 177001(2009).
J. H. Shim, K. Haule, and G. Kotliar, Phys. Rev. B ,060501(R) (2009). J. K. Wang, L. L. Zhao, Q. Yin, G. Kotliar, M. S. Kim,M. C. Aronson, and E. Morosan, Phys. Rev. B , 064428(2011). J. Park, G. Lee, F. Wolff-Fabris, Y. Y. Koh, E. J. Eom,Y. K. Kim, M. A. Farhan, Y. J. Jo, C. Kim, J. H. Shim,and J. S. Kim, Phys. Rev. Lett. , 126402 (2011).
G. Lee, M. A. Farhan, J. S. Kim, and J. H. Shim, Phys.Rev. B , 245104 (2013). Y. J. Jo, J. Park, G. Lee, M. J. Eom, E. S. Choi, J. H.Shim, W. Kang, and J. S. Kim, Phys. Rev. Lett. ,156602 (2014).
L.-L. Jia, Z.-H. Liu, Y.-P. Cai, T. Qian, X.-P. Wang,H. Miao, P. Richard, Y.-G. Zhao, Y. Li, D.-M. Wang,J.-B. He, M. Shi, G.-F. Chen, H. Ding, and S.-C. Wang,Phys. Rev. B , 035133 (2014). Y. Feng, Z. Wang, C. Chen, Y. Shi, Z. Xie, H. Yi,A. Liang, S. He, J. He, Y. Peng, X. Liu, Y. Liu, L. Zhao,G. Liu, X. Dong, J. Zhang, C. Chen, Z. Xu, X. Dai,Z. Fang, and X. J. Zhou, Sci. Rep. , 5385 (2014). K. Wang and C. Petrovic, Phys. Rev. B , 155213 (2012). K. Wang, D. Graf, and C. Petrovic, Phys. Rev. B ,235101 (2013). X. Shi, P. Richard, K. Wang, M. Liu, C. E. Matt, N. Xu,R. S. Dhaka, Z. Ristic, T. Qian, Y.-F. Yang, C. Petrovic,M. Shi, and H. Ding, Phys. Rev. B , 081105(R) (2016). N. Katayama, K. Kudo, S. Onari, T. Mizukami, K. Sug-awara, Y. Kitahama, K. Iba, K. Fujimura, N. Nishimoto,M. Nohara, and H. Sawa, J. Phys. Soc. Jpn. , 123702(2013). H. Yakita, H. Ogino, T. Okada, A. Yamamoto, K. Kishito,T. Tohei, Y. Ikuhara, Y. Gotoh, H. Fujihisa, K. Kataoka,H. Eisaki, and J. ichi Shimoyama, J. Am. Chem. Soc. , 846 (2014).
C.-J. Kang, T. Birol, and G. Kotliar, Phys. Rev. B ,014511 (2017). J. W. Harter, H. Chu, S. Jiang, N. Ni, and D. Hsieh,Phys. Rev. B , 104506 (2016). Y. Mizuguchi, Y. Hara, K. Deguchi, S. Tsuda, T. Yam-aguchi, K. Takeda, H. Kotegawa, H. Tou, and Y. Takano,Supercond. Sci. Technol. , 054013 (2010). X. Wu, C. Le, Y. Liang, S. Qin, H. Fan, and J. Hu, Phys.Rev. B , 205102 (2014). X. Wu, S. Qin, Y. Liang, C. Le, H. Fan, and J. Hu, Phys.Rev. B , 081111(R) (2015). Y. Yasui, H. Sasaki, M. Sato, M. Ohashi, Y. Sekine,C. Murayama, and N. Mˆori, J. Phys. Soc. Jpn. , 1313(1999). E. J. T. Salter, J. N. Blandy, and S. J. Clarke, Inorg.Chem. , 1697 (2016). V. M. Zainullina and M. a. Korotin, Phys. Solid State ,978 (2011). V. M. Zainullina, N. A. Skorikov, and M. A. Korotin,Phys. Solid State , 1864 (2012). P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kavasnicka,and J. Luitz,
Wien2k (Karlheinz Schwarz, Tech. Univer-sit¨at Wien, Vienna, 2001).
S. Qin, Y. Li, Q. Zhang, C. Le, and J. Hu, Frontiers inPhysics (2018), https://doi.org/10.1007/s11467-018-0745-7, arXiv:1702.08304. H. He and M. Aaronson, “private communication,”(2017).
C.-J. Kang and G. Kotliar, Phys. Rev. Materials , 034604(2018). M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev.Lett. , 14 (2007). A. L. Kutepov, Phys. Rev. B , 195120 (2017). C. Franchini, A. Sanna, M. Marsman, and G. Kresse,Phys. Rev. B , 085213 (2010). W. Kang and M. S. Hybertsen, Phys. Rev. B , 195108(2010). C. H. P. Wen, H. C. Xu, Q. Yao, R. Peng, X. H. Niu,Q. Y. Chen, Z. T. Liu, D. W. Shen, Q. Song, X. Lou,Y. F. Fang, X. S. Liu, Y. H. Song, Y. J. Jiao, T. F.Duan, H. H. Wen, P. Dudin, G. Kotliar, Z. P. Yin, andD. L. Feng, Physical Review Letters (2018),https://doi.org/10.1103/PhysRevLett.121.117002,arXiv:1802.10507.
B. Fultz, Prog. Mater. Sci. , 247 (2010). K. Haule and T. Birol, Phys. Rev. Lett. , 256402(2015).
K. Haule and G. L. Pascut, Physical Review B , 195146(2016), arXiv:1602.02819. X. Deng, L. Wang, X. Dai, and Z. Fang, Phys. Rev. B , 075114 (2009). Y. X. Yao, J. Schmalian, C. Z. Wang, K. M. Ho, andG. Kotliar, Phys. Rev. B , 245112 (2011). N. Lanata, X. Deng, and G. Kotliar, Phys. Rev. B ,081108(R) (2015). W.-S. Wang, X.-M. He, D. Wang, Q.-H. Wang, Z. D.Wang, and F. C. Zhang, Phys. Rev. B , 125105 (2010). M. Sandri, M. Capone, and M. Fabrizio, Phys. Rev. B , 205108 (2013). P. Liu, M. Kaltak, J. Klimeˇs, and G. Kresse, Phys. Rev.B , 165109 (2016). H. Peng and S. Lany, Phys. Rev. B , 174113 (2013). J. A. Schiller, L. K. Wagner, and E. Ertekin, Phys. Rev.B , 235209 (2015). T.-H. Lee, Y.-X. Yao, V. Stevanovi´c, V. Dobrosavljevi´c,and N. Lanat`a, arXiv:1710.08586 (2017).
E. B. Isaacs and C. Wolverton, Phys. Rev. Materials ,063801 (2018). J. Harl and G. Kresse, Phys. Rev. Lett. , 056401(2009).
D. G. Schlom, L.-Q. Chen, X. Pan, A. Schmehl, and M. A.Zurbuchen, Journal of the American Ceramic Society ,2429 (2008). “The nomad repository http://nomad-repository.eu/,” . Y. Li, J. Hao, H. Liu, Y. Li, and Y. Ma, J. Chem. Phys. , 174712 (2014).
S. Lany, Phys. Rev. B , 245207 (2008). L. Wang, T. Maxisch, and G. Ceder, Phys. Rev. B ,195107 (2006). A. Jain, G. Hautier, S. P. Ong, C. J. Moore, C. C. Fischer,K. A. Persson, and G. Ceder, Phys. Rev. B84