Quantum Monte Carlo analysis of a charge ordered insulating antiferromagnet: the Ti 4 O 7 Magnéli phase
Anouar Benali, Luke Shulenburger, Jaron T. Krogel, Xiaoliang Zhong, Paul R. C. Kent, Olle Heinonen
QQuantum Monte Carlo analysis of a charge orderedinsulating antiferromagnet: the Ti O Magnéli phase
Anouar Benali, † Luke Shulenburger, ‡ Jaron T. Krogel, ¶ Xiaoliang Zhong, § Paul R.C. Kent, (cid:107) and Olle Heinonen ∗ , § , ⊥ Argonne Leadership Computing Facility, Argonne National Laboratory, Argonne, Illinois 60439USA, HEDP Theory Department, Sandia National Laboratories, Albuquerque, New Mexico87185 USA, Materials Science and Technology Division, Oak Ridge National Laboratory, OakRidge, Tenessee 37831 USA, Material Science Division, Argonne National Laboratory, Argonne,Illinois 60439 USA, Center for Nanophase Materials Sciences and Computer Science andMathematics Division, Oak Ridge National Laboratory, Oak Ridge, Tenessee 37831 USA, andNorthwestern-Argonne Institute for Science and Engineering, Northwestern University, 2145Sheridan Rd., Evanston, Illinois 60208, USA
E-mail: [email protected]
Abstract
The Magnéli phase Ti O is an important transition metal oxide with a wide range of ap-plications because of its interplay between charge, spin, and lattice degrees of freedom. At ∗ To whom correspondence should be addressed † Argonne Leadership Computing Facility, Argonne National Laboratory, Argonne, Illinois 60439 USA ‡ HEDP Theory Department, Sandia National Laboratories, Albuquerque, New Mexico 87185 USA ¶ Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tenessee 37831 USA § Material Science Division, Argonne National Laboratory, Argonne, Illinois 60439 USA (cid:107)
Center for Nanophase Materials Sciences and Computer Science and Mathematics Division, Oak Ridge NationalLaboratory, Oak Ridge, Tenessee 37831 USA ⊥ Northwestern-Argonne Institute for Science and Engineering, Northwestern University, 2145 Sheridan Rd.,Evanston, Illinois 60208, USA a r X i v : . [ c ond - m a t . m t r l - s c i ] J un ow temperatures, it has non-trivial magnetic states very close in energy, driven by electronicexchange and correlation interactions. We have examined three low-lying states, one ferro-magnetic and two antiferromagnetic, and calculated their energies as well as Ti spin momentdistributions using highly accurate Quantum Monte Carlo methods. We compare our resultsto those obtained from density functional theory-based methods that include approximate cor-rections for exchange and correlation. Our results confirm the nature of the states and theirordering in energy, as compared with density-functional theory methods. However, the energydifferences and spin distributions differ. A detailed analysis suggests that non-local exchange-correlation functionals, in addition to other approximations such as LDA+U to account forcorrelations, are needed to simultaneously obtain better estimates for spin moments, distribu-tions, energy differences and energy gaps. Introduction
Electronic correlations play an important role in compounds containing elements with partiallyfilled d- and f-shells, such as transition metals in rows three and four of the periodic table. Corre-lation effects are responsible for a number of phenomena of intense interest during the past severaldecades.
In particular, the oxides of 3d transition metals, such as Mn, Ti, Fe, Co, Ni, and Zn,exhibit strong electronic correlations because of their rather well localized d-orbitals with manydifferent oxidation states. As a consequence, the 3d transition metal oxides exhibit a wide rangeof phenomena, such as charge-transfer or Mott insulating antiferromagnets, metal-insulator tran-sitions, transitions between paramagnetic and non-trivial magnetic phases, and charge and bonddisproportionation, that originate in coupling between charge, spin and lattice degrees of free-dom.
This has led to intense research over the past decade to create functional materials inwhich transition metal oxides are cornerstones. How to efficiently and accurately capture corre-lation effects in computational modeling of materials systems has long been a crucial but elusiveproblem. This is important not only from a fundamental point of view, but also necessary in or-der to create and control functional materials in which electronic correlations play a fundamental2ole. The electronic correlations in 3d transition metal compounds are often strong enough thatpredictive modeling using standard tools such as Density Functional Theory (DFT) can givequalitatively and quantitatively wrong results. Titanium oxides are one class of versatile transitionmetal compounds that are of interest for a number of applications ranging from resistive memoriesto photocatalysis and dyes in paints.
Because of the many oxidation states of Ti, there exista whole series of Ti oxides, starting with TiO , which is an insulator with a band of about 3 eVwith a Ti-charge of +
4, obtained by successively reducing the average charge of the Ti cations withconcomitant removal of oxygen. An interesting subset of Ti oxides consists of the ordered Magnéliphases Ti n O ( n − ) . This family of compounds exhibits sharp metal-insulator transitions associatedwith pronounced charge and/or orbital ordering. In particular, Ti O is a seemingly innocuouscompound with applications in resistive switching, oxide-based fuel cell electrodes, andin catalysis. Ti O exhibits a complex behavior driven by electronic correlations: It undergoesa phase transition between a non-magnetic conductor and semiconductor upon lowering tempera-ture to 150 K, and a upon further reduction of temperature a transition to a magnetic semiconductorat 120 K, and it has a surprisingly rich energy surface in its low-temperature phase, with com-peting ferromagnetic and antiferromagnetic states. Experimental studies by Marezio et al. using single-crystals in X-ray investigation on Ti O at the transition from semiconductor-metalto semiconductor-semiconductor, and by Lakkis et al. combining X-ray diffraction with elec-tron paramagnetic-resonance studies, arrived at two competing models describing the nature of theelectronic structure of Ti O . The first suggests there must be long-range order for the Ti + ionsand, consequently, ordered bipolarons. The second one suggested that all the bipolarons, whichwere well ordered in the low-temperature phase, remained in the intermediate phase, but becomedisordered. At high temperature, both models suggested the metallic behavior was due to a delo-calization of the 3d electrons.Because of its simple composition, being a binary transition metal oxide, and wide versatility,Ti O is an important and representative compound serving as a testbed to elucidate the role of3lectronic correlations. However, despite the large number of experimental studies of the Ti O -phases at low, intermediate, and high temperatures it has been difficult to theoretically mapout these phases because of the role of the electron-electron interactions beyond mean-field ap-proaches such as the Local (Spin) Density Approximation (LDA) within DFT. Liborio and Har-rison used ab initio thermodynamics, with electronic structure calculations within the DFT+LDAscheme, at low temperatures to predict phase stability in good agreement with experiments and alsoused their model to explain the formation process of the Magnéli phase. Leonov et al. performedan LDA+U calculation on the low temperature phase of Ti O and found an energy gap in goodagreement with experiment and an antiferromagnetic (AF) charge-ordered ground state; Eyert etal. performed LDA band structure calculation on the same low temperature structure but wereunable to reproduce the spin-singlet ground state. More recently, Liborio et al. have performedDFT calculations using DFT with hybrid exchange-correlation (XC) functionals (B3LYP) on low,intermediate and high temperature Ti O suggesting an antiferromagnetic (AF) charge-orderedsemiconducting state at low temperature. Weissmann and Weht also used LDA+U and found thesame AF ground state as Leonov et al. and Liborio et al. In a more recent study, X. Zhong et al. suggested two energetically competing low-energyanti-ferromagnetic states at low temperature using an atomic self-interaction correction (SIC)scheme within DFT implemented in the SIESTA electronic structure code, as well as more ex-pensive hybrid functional (HSE06) all-electron calculations. One state was the previouslystudied AF state by Liborio et al. , while the other one was a new AF state with lower energy.The energy difference between these two AF states was only 4 meV per formula unit, which issmall enough that one can reasonably question the validity of the conclusions.We will here present calculations of the magnetically ordered low-temperature states of theMagnéli-phase oxide Ti O using highly accurate Quantum Monte Carlo (QMC) simulations witha goal of conclusively addressing the energy ordering of them, as well as the detailed nature of thespin ordering of the magnetic low-temperature states. We will show that QMC can be used to obtain4hese states and their energies, and our results also demonstrate more broadly the application ofQMC to systems with complex energy surfaces that include competing non-trivial magnetic statesdriven by exchange and correlation. Our QMC results confirm the results of Zhong et al. usingLDA+SIC (and also hybrid functional calculations) regarding the ground state and adjacent low-energy states of the low-temperature phase of Ti O . This in and of itself is useful and interesting:we can validate these results obtained from the much less expensive LDA+SIC calculations. Wewill also analyze the results obtained from the QMC calculations and the LDA+U and LDA+SICcalculations, as well as total energy from HSE06 hybrid functional calculations, in order to obtainsome insight into where DFT calculations fail, as well as into what QMC does differently fromLDA+U and LDA+SIC for this compound. The insight gained from this analysis will help improveDFT-based methods and guide their usage.It has long been known that the workhorse of materials modeling, Density Functional The-ory, in its standard implementations, such as the Local Density Approximation or the GeneralizedGradient Approximation (GGA), fails miserably for systems where electronic correlations areimportant. For example, DFT-LDA or DFT-GGA fail to predict the correct ground state for theseemingly simple binary oxide NiO, apart from the consistently under-estimated band gap in in-sulators. Various extensions beyond the LDA or GGA have been devised in order to better accountfor electronic correlations, especially those that arise from more localized electronic orbitals, suchas 3d or 4f orbitals. The LDA+U (or GGA+U, depending on the choice of exchange-correlationfunctional) method assigns a Hubbard-like on-site repulsion for doubly occupied orbitals and hasbeen used rather extensively with some success to model a large range of compounds. The U -parameter leads to split upper and lower Hubbard bands, with a splitting that scales as U , whichis therefore some measure of the strenght of correlations, as represented by the on-site repulsion.Nevertheless, the LDA+U method suffers from being somewhat ad hoc in that it cannot bederived from more fundamental principles. Furthermore, the U -parameter is typically adjusted tofit something, typically a direct band gap, although more recent methods use self-consistencythrough linear response to arrive at a value for U . In any case, there is no guarantee that the5btained band structure is correct, other than the band gap if that was fitted. Another methodconsists of removing self-interactions, as approximate DFT functionals includes interactions ofeach electron with itself. For itinerant electrons, self-interactions do not lead to large errors, butfor localized electronic states self-interactions can lead to significant errors. Based on work byPerdew and Levy, Spaldin and Filipetti developed a SIC scheme using a unitary transforma-tion of localized orbitals to atomic orbitals and subsequent minimization of the self-interactionerror. This scheme has been implemented in SIESTA, and has shown to give good results fora number of transition metal oxides as far as the nature of the ground state and the direct bandgap are concerned. However, this SIC scheme, while having a more ab initio underpinningthan LDA+U, still fundamentally approximates correlation energies. Furhermore, it includes anadjustable parameter α in the range of [ , ] , with α = α = α ≈ . e.g., the band gap).There are other methods more explicitly devised to deal with electronic correlations, but theseare universally much more computationally expensive than DFT methods. The dynamical meanfield theory (DMFT) has been used rather extensively recently to address questions in a varietyof correlated materials, such as the “Fermi arc” in high-T c materials. The dynamical correlationsin DMFT are treated within a quantum point defect (Anderson impurity) approximation, in whichthe self energies depend on frequency but not wavevector, although more recent developments haveextended DMFT to include some spatial variation, and therefore wavevector dependence. The Density Matrix Renormalization Group (DMRG) is in principle an exact method forobtaining the ground state of a quantum system. However, it is most suited for one-dimensionalor pseudo-one-dimensional systems such as quantum spin chains or ladders. It rapidly becomesintractable as the width of the ladders increase, and it is at the present not at all applicable to real6hree-dimensional materials.QMC simulations are also in principle exact and include all electron-electron interactions. TheQMC method is based on a strict variational principle and so obtains a rigorous upper bound forthe ground state energy by directly minimizing the expectation value of the Hamiltonian for amany-body wavefunction. The one fundamental approximation in QMC for fermionic systems isrelated to the antisymmetry of the wavefunction. The most common way to enforce antisymme-try is through the fixed-node approximation, in which the nodes of the variational wavefunctionare fixed based on an initial antisymmetric trial wavefunction. Note that even with this approxi-mation, QMC still obeys the variational principle and yields an upper bound to the exact groundstate energy, but the quality the upper bound depends on the nodal structure of the initial trialwavefunction. QMC is usually implemented in two steps. The first step is an initial variationalMonte Carlo minimization in which Jastrow-type (or dynamical) correlation factors are added toan initial Slater determinant trial wavefunction, and these correlation factors contain some ad-justable parameters over which the energy is minimized. The second step is a diffusion MonteCarlo (DMC) minimization, in which the many-body Schrödinger equation is written as a dif-fusion equation in imaginary time. The expectation value of the Hamiltonian is then minimizedby propagating a convolution of the wavefunction with the Green’s function in imaginary timeusing many random walkers. Due to its numerical expense, QMC methods were for a long timelimited to model systems or small atoms or molecules. However, because the random walkers inthe DMC minimization are only weakly dependent, these methods scale linearly with nearly anynumber of cores available (so long as the memory per core is sufficient). The arrival of leadership-class computers such as the Blue Gene/Q or Cray Titan XK7 with hundreds of thousands of coreshave led to QMC methods being routinely applicable to a wide variety of systems, from simplemolecules to condensed phase metals. More recently, Shulenburger and co-workers applied QMC to study inter-layer interactions in bulk and few-layer phosphorous and showed thatDFT-based methods including dispersive forces cannot account for a charge redistribution thatarises because of inter-layer interactions. Here, we use QMC to provide an upper bound and error7stimates for the LT-phases of Ti O , with 22 atoms in the unit cell and several different magneticstates, and to analyze the differences between the essentially exact QMC results and those obtainedfrom LDA+U and LDA+SIC. As a side result, we also show how DMC can be used in an un-biased way to select a value for the Hubbard U -parameter in an LDA+U scheme. Methods
QMC calculations
This study was performed using DMC within the fixed node approximation. Since QMC, andmore specifically, DMC, have been thoroughly reported, we will only briefly describe themain features of our implementation in the QMCPACK package.
To describe our systems we used a single-determinant Slater-Jastrow wave function Ψ ( x , ..., x N ) expressed as Ψ ( x , ..., x N ) = e J ( x ,..., x N ) Ψ AS ( x , ..., x N ) (1)where x i ≡ { r i , σ i } is a space-spin coordinate, J ( x , ..., x N ) is the Jastrow factor describing the cor-relation between electrons, Ψ AS ( x , ..., x N ) is the antisymmetric fermionic part of the wavefunctionand N is the total number of electrons.The antisymmetric part is calculated within the framework of the DFT, using the plane wavepseudopotential code PWSCF contained in the Quantum ESPRESSO (QE) simulation package. We used the rotationally invariant LDA+U implementation in QE with the LDA XC functionalproposed by Perdew and Zunger and a single Slater coefficient F for the U -parameter. This isused as an additional variational optimization parameter. Convergence with respect to the kineticenergy cutoff ( E w f ccut ), the density cutoff ( E rhocut ) and size of Monkhorst-Pack mesh of k -pointswere checked in order to achieve the same precision in all calculations, leading to the followingvalues: E w f ccut =
300 Ry, E rhocut = × × k -point grid in the HCP primitive8ell for the Ti bulk calculations. The single-atom calculation calculations were performed in cubicboxes with 30 a.u dimension ensuring an insignificant amount of interaction between atoms andtheir periodic images.The Jastrow factor can be factored to one-body (effective electron-ion interactions), two-body(electron-electron interactions), three-body (electron-electron-ion interactions) and higher-bodyorders ( J = J J J .. ) . In this work, using one-body and two-body terms provided a satisfactorytrade-off between computational expense and ease of optimization: J ( r , ..., r N ) = ∑ iI χ I ( r iI ) − ∑ i < j u ( r ij ) , (2)where χ and u are one-dimensional optimized splines; i labels the electron position, I labels theion position, r i j = | r i − r j | is the electron-electron distance and r iI = | r i − R I | is the electron-iondistance. In this definition, χ is the cusp-less electron-Ion Jastrow factor, while u describes theelectron-electron correlations. Different correlation factors were used for like and unlike spins.All the Jastrow parameters were optimized within the energy minimization framework by min-imizing a combination of the total energy and variance. Since the Jastrow prefactors are strictlypositive, the nodes of the Slater-Jastrow wave function are determined by ψ AS and therefore thefixed-node error in the DMC calculations is entirely due to the Slater determinant obtained fromthe DFT calculation.All calculations were performed within the pseudopotential (PP) approximation, in which thecore electrons are replaced by a nonlocal potential operator. In the DMC part, we used the lo-cality approximation to evaluate the nonlocal PP. Norm-conserving pseudopotentials (NCPP)were used for both QMC and DFT calculations. They were defined with ( s p d s ) valenceelectrons for Titanium atoms and ( s p ) for Oxygen atoms and generated with the OPIUM code. The oxygen PP was generated and optimized by Shulenburger and Mattsson. The Ti PPwas generated by Krogel et al. using the Opium code following Rappe et al. scheme. The9any-body wavefunction was evaluated using real space cubic b-splines for speed, and we useda minimum of 8192 walkers with a 0.005 a.u time step for the DMC calculations with extrapola-tions to the zero timestep limit. Finite-size effects - Solids
When determining properties of solids in the thermodynamic limit with QMC using periodicboundary conditions, it is necessary to control finite-size errors that arise from the spurious cor-relations created by the periodic images of the system. These errors can be grouped into twocategories: single-particle effects and many-body effects. In single-particle theories (and effectivesingle-particle theories such as DFT), calculations are usually reduced to the primitive unit cell.Then the properties of a system in the thermodynamic limit may be determined by using periodicboundary conditions and integrating over the first Brillouin zone. Although such an integration isalso used in a many-body theory such as DMC, there is an additional finite-size error due to imagecell correlations, and more sophisticated techniques are needed to extrapolate to the thermody-namic limit. For this reason, DMC calculations of periodic solids are generally performed usingsupercells to minimize this error. When simulating systems with periodic boundary conditions,the energy contributions from the classical Coulomb interactions cannot be obtained by the sumsover the infinite number of images do not converge when using the classical 1 / r Coulomb potential.The common solution is to use the Ewald summation techniques. However, the periodicity of thesystem implies the periodicity of the electronic correlations, which introduces effects of electronsinteracting with their periodic images in neighboring cells. This unphysical error corresponds to anXC energy of a system with a periodic XC hole. These errors are not present in DFT simulationsbecause the XC interactions are derived from the electron gas converged at the limit of infinitelylarge simulation cells. The first option for finite-size correction consists of increasing the systemsize in order to reduce the spurious effective interactions between an electron and its XC hole. Thisleads to a rapidly increasing cost of the calculations with the number of electrons in the simulationcell. The second finite-size correction is to add the model periodic Coulomb (MPC) interaction to10he Ewald energy as suggested by Fraser et al . A final correction is applied to the kinetic energyfollowing the scheme of Chiesa et al . For all calculations in this study, energy convergence wasreached at simulation cells containing 88 atoms and applying both MPC and Chiesa corrections.In periodic boundary conditions, the Hamiltonian is invariant with respect to translating anyparticle around the periodic boundaries. According to Blochâ ˘A ´Zs theorem, this implies that thatany solution can be characterized by a given vector twist angle θ = ( θ x , θ y , θ z ) , defined by Ψ ( r + L x ˆ x , r , . . . , r N ) = e i θ x Ψ ( r , r , . . . , r N ) , (3)with similar definitions for θ y and θ z . In order to obtain results in the thermodynamic limit, all k -points in the Brillouin zone have to be averaged over. Because a change in the twists changes thegrid over which the kinetic energy is calculated, an average over twist angles is an efficient wayto approach the thermodynamic limit. Specifically, an average is taken of expectation values overall the k -points in the primitive cell. DMC calculations at different twist angles are statisticallyindependent, so the cost of twist averaging is negligible except for initialization and equilibration.The error can therefore be reduced by increasing the density of the twist grid. To take into accountthe two-body-finite-size effect, we increased the size of the supercell until the change in energyfrom one supercell size to the next is equal to 0.1mHa/unit-cell, leading to a supercell 4 timeslarger than the primitive cell (Figure-1) and 32 boundary twist per supercell.Ti O at low temperature (less than 120 K) was modeled using two formula unit cells (Ti O )to account for the various spin configurations, as shown in Fig. 1. The enlarged unit cell thus con-tains one chemical unit cell plus one image of it. In order to distinguish between four Ti atoms andtheir images in the enlarged unit cell, we adopt the following notation for the remaining of the dis-cussion; (Ti1,Ti2, Ti3, Ti4) and their images (Ti1’, Ti2’, Ti3’, Ti4’). The atomic spins around theTitanium atoms Ti2, Ti4, and their periodic images Ti2’, Ti4’, in the simulation cell are very small.Zhong et al , using LDA+SIC as well as the HSE06 hybrid functional, identified a ferromagneticstate (FM), where all Ti atom have a spin up ( ↑ ), and three anti-ferromagnetic state, AF1 AF2 and11igure 1: Simulation cell of (Ti O ) containing 22 atoms (14 oxygen and 8 titanium)Table 1: Spin configuration for Titanium atoms in Ti O ↑ ) represents a spin up while ( ↓ ) represents a spindown. The spins are for all Ti x atoms and their images Tix’ in the sublattice as follow (Ti1 , Ti2,Ti3, Ti4 | Ti1’ ,Ti2’ ,Ti3’ ,Ti4’); 0 denotes an almost null atomic spin.Spin ConfigurationFM ( ↑ , ↑ , ↑ , ↑ | ↑ , ↑ , ↑ , ↑ )AF1 ( ↑ , , ↓ , | ↑ , , ↓ , ↑ , , ↑ , | ↓ , , ↓ , ↑ , , ↓ , | ↓ , , ↑ , Results and Discussion
As we discussed in previous sections, the results obtained using DFT are highly dependent onthe choice of the XC functionals. In the case of the antiferromagnetic states of Ti O , none of12he tested functionals (LDA or GGA) was able to reproduce or predict the AF states (they werenot obtained as stable states). Using LDA+U, LDA+SIC, or the HSE06 hybrid functional, theAF states are stable but the difference in energy is only 4 - 5 meV per formula unit. In previousstudies using LDA+U, the value of U was set to reproduce the experimental band gap of Ti O ( U = . J = . used a value of U = .
44 eV obtained from considering energy differences between phases and obtained a bandgap of 0.7 eV. In our LDA+U calculations, we used LDA+U as implemented in the QuantumEspresso code with a single Slater coefficient F = U . To select the value of U in our study, wechose to use the variational nature of QMC to scan through the values U . Doing so, we remove thead-hoc empirical contribution to the Hamiltonian. As expected and shown on Fig. 2(a), increasingthe value of U from 1 to 10, which is equivalent to populating the d band, leads to a linear increaseof the energy at the DFT level. In contrast, the total DMC energy (as well as the VMC one) ofthe AF1 state as a function of the trial wavefunction, characterized by the U -parameter, follows anonlinear U dependence with a minimum at U DMC = . ± . U -dependent trial wavefunction.Thisshows clearly the sensitivity of the energetics to the U parameter as more sophisticated treatmentsof dynamic correlation are used; whereas the variation in the value of U between 1 eV and 10 eVleads to a spread of 20 eV in LDA+U, it only leads to a spread of 1.5 eV in DMC. In the remainingof the calculation we used the U DMC = . It is interesting to note that thisunique method selection of U is quite different from more recent LDA+U methodologies, where U is selected within a rotationally invariant simplified LDA+U scheme (keeping only the F Slatercoefficient) by enforcing self-consistency between the chosen U and the value of U obtained fromlinear response. Within that scheme, the total energy is still a monotonic function of U , in contrastwith the dependence of the energy on U obtained here from DMC. It should be pointed out thatthere does not exist, to the best of our knowledge, any fundamental justification for either methodof selecting U . 13 ohesive Energy evolution of Ti4O7 at 120K in the AF1 state with U-105-100-95-90-85-80-75 0 2 4 6 8 10 E ne r g y ( e V ) (a) Value of U (eV) DFT-103.5-103-102.5-102-101.5-101-100.5-100-99.5 0 2 4 6 8 10 E ne r g y ( e V ) (b) Value of U (eV)Minimum E (DMC) = -103.003 +/- 0.039 eV Minimum U (DMC) = 4.327 +/- 0.297 eV DMC Figure 2: DFT (a)and DMC (b) energies of Ti O as a function of the value of U for the AF1state. The DFT curves increases linearly with the value of U while DMC one is non linear functionshowing a minimum in energy at U = . using α = .
5, and also from HSE06 hy-brid functional; Table 3 shows the energy gaps obtained from LDA+U, LDA+SIC calculations( α = . U -parameter in LDA+U or the α -parameter in LDA+SIC in our calculations was adjusted to fit the energy gaps. The LDA+U,LDA+SIC, and HSE06 hybrid functional calculations show a clear separation in energies betweenthe FM state and the AF states, especially the LDA+U and HSE06 calculations, while the energydifference between AF1 and AF3 is very small (4 - 5 meV), clearly beyond errors within LDA+U,and LDA+SIC. However, the DMC results increase significantly the energy differences betweenthe AF states, favoring the AF3 state (staggered antiferromagnetic order within each sublattice andbetween sublattices) over both the AF1 (staggered antiferromagnetic order within each sublattice,but ferromagnetic order between sublattices) and the FM states, although the DMC energy sepa-ration is substantially smaller between the FM and AF3 state than obtained from LDA+U. Notethat the DMC calculations also have error bars and the energy differences between the states aresignificantly larger than the error bars. It is important to note that this ordering is identical to whatis found with DFT; the main difference is the accuracy of the prediction. The energy difference14ound between FM and AF3 at the DMC level of theory becomes much closer to the LDA+SIC re-sult despite starting from a much larger difference in the LDA+U calculations which superficiallysuggests a somehow a better correction to the XC interactions for LDA+SIC than for LDA+U. Wewill argue later that this is not necessarily the case.We have also calculated the energy differences between the AF3 and AF1 states for a rangeof values of the U -parameter in LDA+U, for the α -parameter in LDA+SIC, and for the exchangemixing parameter a in the HSE hybrid functional. We note that recent work shows that the choicesof the exchange mixing of a = .
303 and range ω = .
201 Å − specified in HSE06 are optimal. The energy difference between AF3 and AF1 is relatively insensitive to the values of U , α , orexchange mixing a . As U , α , or a increase, the 3d states become more localized, and a reducedoverlap leads to a smaller energy difference. At the same time the band gap increases by pushingapart the split bands in LDA+U, or by reducing the energy of the occupied orbitals relative to theunoccupied ones in LDA+SIC, or by increasing the exchange interactions between 3d states as a increases. When U varies from 1 eV to 7 eV, the energy difference from LDA+U only varies from4.1 meV to 1.6 meV. Similarly, when α varies from 0.3 to 0.7, the energy difference varies from4.5 meV to 2.3 meV, and when a varies from 0.15 to 0.4, the energy difference decreases from5.3 eV to 3.3 meV. On the other hand, the band gaps increase approximately linearly from about0.3 eV to about 2.8 eV for LDA+U, from about 0.2 eV to 1.7 eV for LDA+SIC, and from about0.5 eV to 3 eV for the HSE functional. Moreover, as U, α , or a decrease, the local moments de-crease because of delocalization of the orbitals. For example, the (Mulliken population) momentsfrom LDA+SIC decrease in magnitude from 0.93 to 0.68 as α decreases from 0.7 to 0.3, and theorbital moments in HSE decrease from 0.96 to 0.76 as a decreases from 0.4 to 0.15: improving theband gap by decreasing α or a makes the magnetic moments worse while the energy differencebetween AF3 and AF1 barely changes.Figure 3 shows the DMC iso-surfaces of the spin densities with yellow representing net up-spin(positive spin density) and blue representing down-spin (negative spin density). We note first of all15igure 3: Spin density difference ( ρ ↑ − ρ ↓ ) from the DMC (top), LDA+U (center) and LDA+SIC(bottom) depicting iso-surfaces of 0.0298 electrons per a , where a is the Bohr radius, in the FMstate (left panel), the AF3 state (middle panel), and the AF1 state (right panel). The value of theiso-surfaces was chosen to give a good visual representation of the spin densities around the Ti + sites. Yellow represents up-spin ( ρ ↑ − ρ ↓ = . / a ), blue down-spin ( ρ ↑ − ρ ↓ = − . / a )16able 2: Energy difference between the antiferromagnetic state (AF3) and the ferromagnetic state(FM), and the antiferromagnetic state (AF2) and the ferromagnetic state (FM) in eV per formulaunit using LDA+U, LDA+SIC and DMC. A positive value means that the AF3 state is lower inenergy. Method E(FM-AF3) E(AF1-AF3)LDA + U +0.608 +0.004LDA + SIC +0.118 +0.004HSE06 +0.605 +0.005DMC +0.17(1) +0.07(1)Table 3: Energy gaps in eV obtained from LDA+U, and LDA+SIC calculations for the FM, AF1,and AF3 states. The experimental gap is about 0.25 eV.Method FM AF1 AF3LDA + U 1.32 1.55 1.61LDA+U d yz spatial distribution (with the x-axisthe projected line connecting Ti3 - Ti1 - Ti1’- Ti3’) in all three magnetic states – the momentsonly change sign between the different magnetic states. There is only a very slight rotation ofthe moment distributions on Ti3 and Ti3’ in the AF state that make them appear like a linearcombination of d yz and d x − y orbitals. In contrast, the spin distributions on 3 and 3’ have differentshapes in the FM state compared to in the AF states and are neither pure d z or four-lobed planar d -orbitals.For comparison, the middle row of Fig. 3 depicts iso-spin density surfaces obtained fromLDA+U. This figure shows that in LDA+U too the shapes of the spin density distributions inAF3 (center) and AF1 (right) are identical, and extremely similar to those obtained from DMC,with d yz character on Ti1 and Ti1’, and a combination of d yz and d x − y on Ti3 and Ti3’, also ingood qualitative agreement with the results of Leonov et al. In contrast, the spin distributions inthe FM state show some marked differences both from those in the AF states as well as those in theDMC FM state.
All spin distributions in FM have approximately the same d yz + d x − y character,with spin density lobes much more directed along the line connecting the Ti3 - Ti1 - Ti1’ - Ti3’17ites. This clearly leads to a different near-neighbor overlap in the LDA+U FM than in the DMCFM state.The lower panels of Fig. 3 show iso-surfaces of the spin densities obtained from the LDA+SICcalculations. The spin distributions in AF3 and AF1 on Ti1 and Ti1’ again agree very well with theDMC and LDA+U results. The shapes of the distributions at Ti1 and Ti1’ and also at Ti3 and Ti3’are completely different in the FM state between DMC and LDA+SIC. Now the Ti1 and Ti1’ siteshave much more d xz character, while the Ti3 and Ti3’ sites have d yz + d x − y character. Again, thedifferences in spin distributions obtained from LDA+SIC compared to DMC will obviously leadto different overlaps between near neighbor sites.Before we discuss these results in detail, it is worth recapitulating the basic ingredients inLDA+U and LDA+SIC. First, both these methods as usually implemented involve projecting thesingle-particle states on a local atomic basis. In LDA+U, a Hubbard-like energy that penalizes twoelectrons (spin-up and spin-down) on the same atomic orbital is added to the total energy. Anotherterm that (approximately) corrects for double-counting also has to be added to the total energy.The effective potential added to the Kohn-Sham equations then contains a spin-independent anda spin-dependent part, in addition to the standard LDA or GGA potentials. Note that this con-struct projects out self-interactions between the different angular momentum atomic orbitals. Thespin-dependent part, together with the LDA or GGA potentials, give rise to an effective exchangeinteraction between neighboring sites. This interaction comes from a local spin-dependent po-tential that arises from the LDA approximation and the on-site U -interaction. Therefore, the near-neighbor exchange is a consequence of the single-particle states not being perfectly localized oneach atomic site, but having finite overlaps with atomic orbitals on neighboring sites. Nevertheless,LDA+U has been successful in obtaining a reasonable ground state for a number of antiferromag-netic Mott insulators or charge transfer insulators, in addition to obtaining reasonably accuratevalues for spin moments and effective exchange couplings. The SIC implementation also projects the single-particle states onto localized orbitals byseeking a unitary transformation to localized orbitals that minimizes the self-interaction energy18etween these orbitals. Main differences between SIC and LDA+U are that SIC only operates onoccupied orbitals and pushes them down in energy, while the unoccupied orbitals are not affected,and that SIC does not add additional ad hoc interactions to the Hamiltonian; SIC only tries toeliminate self-interactions. Both LDA+U and SIC correct the tendency of local or semi-local XCfunctionals to erroneously delocalize d-orbitals of transition metals. This delocalization, which islargely due to incorrect self-interactions of the orbitals, leads to overestimating the Heisenberg cou-pling between local moments. Interestingly, the self-interaction free Hartree-Fock approximationtends to over-correct and leads to too small effective Heisenberg couplings, presumably because ofa lack of any correlation energy. An effective Heisenberg interaction arises in SIC, as in LDA+U,from the overlap of atomic orbitals on neighboring sites, but the strength of the interaction nowcomes only from projecting out self-interactions (in addition to the local LDA xc-contribution)and has no contribution from an added local interaction, as in LDA+U. The SIC implementation inSIESTA has also been used to calculate spin moments as well as effective Heisenberg exchange ina number of structures, but the results for the Heisenberg couplings are mixed with relatively goodresults for structures with half-filled d-bands but worse away from half filling. Notably, awayfrom half filling, SIC appears to need different values of the α -parameter depending on whichquantity is sought, the band gap or the exchange couplings.With this background we now turn to discussing the results in more detail. We first calculate thelocal spin moments associated with each Ti atom. The idea is that the energy differences betweenthe states can be analyzed using a Heisenberg model, and this analysis can therefore also givesome insight into how LDA+U and LDA+SIC treat electron-electron interactions differently andincorrectly. Figures 4 and 5 shows the calculated spin moments obtained by integrating the spindensity in a sphere of radius ranging of 1.4 Å centered on the Ti atoms in the unit cell. Note thatthe spin moments on Ti1, Ti1’, Ti3, and Ti3’ are very nearly equal in magnitudes in the FM, AF1,and AF3 states, while the spin moments on the remaining Ti ions are close to zero ( ∼ . r = .
802 Å, r = .
133 Å, and r = .
159 Å between Ti1 and Ti3 (and Ti1’ and Ti3’), Ti1 and Ti1’ (and Ti3 and Ti3’), and Ti1’19igure 4: Calculated spin moments from DMC (black squares), LDA+U (red triangles), andLDA-SIC (blue diamonds) for AF3 (left panels) and AF1 (right panels). The bottom row showsthe moments on an enlarged scale. The moments were obtained by integrating spheres of radius1.4 Å centered on the Ti atoms. 20igure 5: Calculated spin moments from DMC (black squares), LDA+U (red triangles), andLDA-SIC (blue diamonds) in the FM state in units of the Bohr magneton. The right panel showsthe moments on an enlarged scale. The moments were obtained by integrating spheres of radius1.4 Å centered on the Ti atoms.and Ti3 (and Ti1 and Ti3’). This lends justification to the hypothesis that the energy differences areto lowest order driven by nearest-neighbor exchange coupling. First of all, the large spin momentson Ti1, Ti3, Ti1’, and Ti3’ are nearly identical in magnitude both for DMC, LDA+U, and LDA-SIC in AF3, in FM, as well as in AF1; the FM moments from LDA+SIC are only slightly largerin the FM state. This is in and of itself rather remarkable: both LDA+U and SIC seem to handlecharge disproportionation and local spin moments very well. LDA+U captures very accurately thelarge spin moments on Ti1, Ti1’, Ti3, and Ti3’ in all three phases. It also does a good job withthe small spin moments on Ti2, Ti2’, Ti4, and Ti4’ in the FM and AF1 states; the main errors arethe Ti4 and Ti4’ moments in the AF3 state, with approximately twice as large magnitude for thespin moments as obtained from DMC. LDA+SIC does a little worse in terms of the spin moments,with the large moments slightly overestimated in the FM state, and the small Ti4 and Ti4’ momentmagnitudes much overestimated in the AF3 state (about 0.03 Bohr magnetons instead of 0.005 -0.01), and slightly overestimated magnitudes of the small moments Ti2, Ti2’, Ti4 and Ti4’ in theAF1 state. For comparison, Weissman and Weht obtained a spin moment of magnitude 0.64 µ B on the large-moment Ti sites, and Leonov et al. obtained 0.67 µ B , and a magnitude of between0.02 and 0.04 µ B for the small-moment Ti sites.21e first focus on the FM and AF1 states that are separated by 0.1 eV in the DMC, comparedto 0.604 eV in LDA+U, and 0.114 eV in LDA+SIC. Because of the identical spin distributions inDMC on Ti1 and Ti1’ in the FM and AF1 states, the energy difference between these states mustbe driven mainly by an effective antiferromagnetic exchange coupling between the pairs nearest-neighbor pairs Ti1-Ti3 and Ti1’-Ti3’. Note now the small but clearly visible difference in the spindistributions on the Ti3 and Ti3’ sites between the FM and AF1 states in both the DMC resultsand the LDA+U (Fig. 3 middle panels). The orientation of the lobes of the spin distributions in theLDA+U FM clearly suggests a larger overlap between nearest neighbors Ti1 - Ti3 and Ti1’ - Ti3’than in DMC and a concomitant overestimated exchange coupling between these sites, consistentwith a much too large energy difference between AF1 and FM in LDA+U, compared to DMC.The picture In LDA+SIC is less clear with the near-perpendicular orientation of the lobes on Ti1and Ti3, and Ti1’ and Ti3’ in FM, compared to DMC but visual inspection of Fig. 3, lower leftpanel, makes it plausible that the overlaps between the Ti1 and Ti3 sites, and the Ti1’ and Ti3’ sitesare slightly larger in LDA+SIC than DMC, with a slightly larger exchange coupling and energydifference FM-AF1.LDA+U severely overestimates the Ti1 - Ti3, and the Ti1’ - Ti3’ exchange coupling, whileLDA+SIC slightly overestimates it. Increasing U more strongly localizes the 3d-like electronicstates on single sites and reduces overlap with adjacent sites, which would reduce the effectiveHeisenberg exchange and therefore reduce the energy difference between the FM and AF3 states(neglecting the change in kinetic energy because of reduced near-neighbor hopping). It would alsoreduce the energy gap, which is significantly overestimated with this value of U . But then theatomic moments (and charges) would be overestimated. This points to a problem, at least withthe implementation of LDA+U used here: the value of U is “right" for the localized moments butgives rise to too strong near-neighbor exchange and a much too large energy gap. Furthermore,the shape and orientations of the spin distributions in the FM state are not right in LDA+U, and noadjusting of the U -parameter can correct that.In contrast, LDA+SIC appears to localize the moments slightly too much in the FM state,22ith moments about 10% too large. This is reminiscent of the general trend of the Hartree-Fockapproximation, where exchange-only leads to too strong near-neighbor same-spin repulsion andorbital localization; some corrections for electronic correlations are needed. However, the energygap from LDA+SIC is also too large (see Table 3). Decreasing α will reduce the localization of the3d-orbitals, which increases the near-neighbor Heisenberg exchange and reduces the energy gap. If we use α = . α = .
5, the energy gap is reduced to about 0.2 eV for the AF1 andAF3 states, but now the FM state is metallic. The energy differences become 9 meV betweenAF1 and AF3, and 128 meV between FM and AF3; the energy differences increase somewhat asthe moments delocalize more with reduced self-interaction corrections. The problem is just thatthe moments have delocalized too much with and almost 20% reduction of the large momentsin the AF1 and AF3 states and a somewhat smaller reduction in the FM state. Also, the shapeand orientations of the spin distributions are not right in the LDA+SIC FM state, and tuning the α -parameter is not going to correct the spin distributions.We note that HSE06 calculations, similarly to LDA+U, gives too large an energy separation be-tween the FM and AF3 states while the energy gap closes in the FM state. This is indicative of tooshort range of the exchange interaction, with the consequence that the orbitals are too delocalizedin the FM state.The differences between the AF3 and AF1 states are much more subtle. First of all, the energydifference is much smaller than between the AF1 and FM states, only about half from DMC, whileLDA+U and LDA+SIC severely underestimates this energy difference at only 1 %-3 % of the AF3-FM one. Within a Heisenberg picture, the energy difference between these two states is driven bynext-nearest neighbor exchange between the pairs Ti1 and Ti1’, and Ti3 and Ti3’. This pictureis supported by the fact that the shapes of the moment distributions are identical (see Fig. 3).However, the large atomic moments are very nearly equal for QMC, LDA+U, and LDA+SIC.Therefore, the energy difference must be driven either by underestimated coupling between thesepairs of moments, or by differences in other pairs of interactions, such as between Ti1 and thesmall moment Ti4, and between Ti1’ and Ti4’. The inter-atomic distance between these pairs of23toms is also 3.133 Å, the same as between Ti1 and Ti1’. Because of this, and because of the verymuch smaller moments in Ti4 and Ti4’, it seems unlikely that the these interactions are relevant:one would expect the overlaps between the Ti1 and Ti1’ orbitals to be approximately the same asbetween the Ti1 and Ti4 atoms so the smaller moments on Ti4 and Ti4’ would make the interactionbetween Ti1 and Ti4 about 1% of the interaction of the exchange interaction between Ti1 and Ti1’.The origin of the difference must therefore come from the interactions between Ti1 and Ti1’, andTi3 and Ti3’. For both LDA+U and LDA+SIC, this implies that either the orbital overlaps aretoo small, or the effective exchange interaction between the orbitals is too small. Given the factthat the spin distributions are very nearly identical in DMC and LDA+U in AF1 and AF3 it seemsmore likely that the LDA+U approximation is too crude and together with the LDA XC functionalsimply does not have the correct form for the XC potential, even if the orbitals are very accurate.This is also consistent with the large energy gaps obtained with LDA+U. In fact, it is hard to seehow a purely local functional, the effective XC functional from the U -parameter and the LDA XCfunctional, can be constructed to simultaneously correct the energy gap and the effective exchangecouplings while keeping the atomic moments, spin distributions, and charges unchanged. Thesame conclusion holds for LDA+SIC, given that the moments and spin distributions in AF1 andAF3 too are very nearly identical to those of DMC. Conclusions
We have here analyzed the low-temperature phase the Magnéli phase Ti O and its non-trivialmagnetic states. We used very accurate Quantum Monte Carlo methods to decisively clarify thenature and energetic ordering of the states. This clearly demonstrates the applicability of highlyaccurate QMC methods to magnetic systems with competing low-energy states. We also comparedour results from DMC with LDA+U and LDA+SIC in order to gain some insight into the errorsgenerated by LDA+U and LDA+SIC. Our results from DMC shows that LDA+U, LDA+SIC, andalso HSE06, give the right states with the right sequence in total energy, while the energy separa-24ions are different from those obtained using DMC; especially the FM energy is overestimated byLDA+U and HSE06 relative to the AF states. The larger energy difference between the AF3 andAF1 states obtained by DMC compared to the results from LDA+U and LDA+SIC is the effectof correctly including electronic correlations. Both LDA+U and LDA+SIC obtain total momentson the Ti + sites in very good agreement with DMC. However, detailed comparisons between theDMC results, on the one hand, and the LDA+U and LDA+SIC results on the other hand, show thatLDA+U and LDA+SIC generate errors in the spin distributions about the Ti + sites in the ferro-magnetic state, while the spin distributions in the AF states are in good agreement with those fromDMC. Furthermore, the energy gap for all three states are severely overestimated both by LDA+Uand LDA+SIC, in spite of the good agreement in the obtained spin moments between the threemethods. This demonstrates that attempts to correct LDA+U and LDA+SIC can in all likelihoodnot be based on local XC functionals: there appear to be no possibilities within a local XC schemeto simultaneously correct the energy gap and the spin spin distributions while at the same timemaintaining the spin moments on the Ti + sites. This suggests that methods that correct LDA+Uor LDA+SIC need to contain more complicated non-local forms of the XC potential. Acknowledgments
An award of computer time was provided by the Innovative and Novel Computational Impact onTheory and Experiment (INCITE) program. This research used resources of the Argonne Leader-ship Computing Facility, which is a DOE Office of Science User Facility supported under ContractDE-AC02-06CH11357. Sandia National Laboratories is a multiprogram laboratory managed andoperated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, forthe U.S. Department of Energy’s National Nuclear Security Administration under Contract No.DE-AC04-94AL85000. AB, LK, JK and PK are supported through Predictive Theory and Model-ing for Materials and Chemical Science program by the Basic Energy Science (BES), Departmentof Energy (DOE). The work by O.H. was supported by the Department of Energy, Office of Sci-25nce, Division of Materials Science and Engineering. X.Z. was supported by U. S. DOE, Office ofScience under Contract No. DE-AC02-06CH11357
References
1. Zaanen, J.; Sawatzky, G. A.; Allen, J. W. Band gaps and electronic structure of transition-metalcompounds.
Phys. Rev. Lett. , , 418–421.2. Anisimov, V. I.; Aryasetiawan, F.; Lichtenstein, A. I. First-principles calculations of the elec-tronic structure and spectra of strongly correlated systems: the LDA + U method. Journal ofPhysics: Condensed Matter , , 767.3. Anisimov, V. I.; Poteryaev, A. I.; Korotin, M. A.; Anokhin, A. O.; Kotliar, G. First-principlescalculations of the electronic structure and spectra of strongly correlated systems: dynamicalmean-field theory. Journal of Physics: Condensed Matter , , 7359.4. Dagotto, E. Complexity in Strongly Correlated Electronic Systems. Science , , 257–262.5. Sarrao, J.; Crabtree, G.; Fleming, G.; Hemminger, J. C.; Ratner, M. Challenges at the Frontiersof Matter and Energy: Transformative Opportunities for Discovery Science. US DOE BasicEnergy Sciences Report ,6. Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas.
Phys. Rev. , , B864–B871.7. Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. , , A1133–A1138.8. Grätzel, M. Perspectives for dye-sensitized nanocrystalline solar cells. Progress in Photo-voltaics: Research and Applications , , 171–185.9. Péchy, P.; Renouard, T.; Zakeeruddin, S. M.; Humphry-Baker, R.; Comte, P.; Liska, P.;Cevey, L.; Costa, E.; Shklover, V.; Spiccia, L. et al. Engineering of Efficient Panchromatic26ensitizers for Nanocrystalline TiO -Based Solar Cells. Journal of the American ChemicalSociety , , 1613–1624, PMID: 11456760.10. Kwon, D.-H.; Kim, K. M.; Jang, J. H.; Jeon, J. M.; Lee, M. H.; Kim, G. H.; Li, X.-S.; Park, G.-S.; Lee, B.; Han, S. et al. Atomic structure of conducting nanofilaments in TiO2 resistiveswitching memory. Nat Nano , , 148–153.11. Hashimoto, K.; Irie, H.; Fujishima, A. TiO 2 Photocatalysis: A Historical Overview and FutureProspects. Japanese Journal of Applied Physics , , 8269.12. Leonov, I.; Yaresko, A. N.; Antonov, V. N.; Schwingenschlögl, U.; Eyert, V.; Anisimov, V. I.Charge order and spin-singlet pair formation in Ti O . Journal of Physics: Condensed Matter , , 10955.13. Yang, Y.; Lu, W. Nanoscale resistive switching devices: mechanisms and modeling. Nanoscale , , 10076–10092.14. Strachan, J. P.; Pickett, M. D.; Yang, J. J.; Aloni, S.; David Kilcoyne, A. L.; Medeiros-Ribeiro, G.; Stanley Williams, R. Direct Identification of the Conducting Channels in a Func-tioning Memristive Device. Advanced Materials , , 3573–3577.15. Tominaka, S. Facile synthesis of nanostructured reduced titanium oxides using borohydridetoward the creation of oxide-based fuel cell electrodes. Chem. Commun. , , 7949–7951.16. Yao, C.; Li, F.; Li, X.; Xia, D. Fiber-like nanostructured Ti O used as durable fuel cell catalystsupport in oxygen reduction catalysis. J. Mater. Chem. , , 16560–16565.17. Wang, L.; Lettenmeier, P.; Golla-Schindler, U.; Gazdzicki, P.; Canas, N. A.; Morawietz, T.;Hiesgen, R.; Hosseiny, S. S.; Gago, A. S.; Friedrich, K. A. Nanostructured Ir-supported onTi O as a cost-effective anode for proton exchange membrane (PEM) electrolyzers. Phys.Chem. Chem. Phys. , , 4487–4495. 278. Bartholomew, R.; Frankl, D. Electrical Properties of Some Titanium Oxides. Phys. Rev. , , 828–833.19. Inglis, A. D.; Page, Y. L.; Strobel, P.; Hurd, C. M. Electrical conductance of crystalline Ti n O n for n=4-9. Journal of Physics C: Solid State Physics , , 317.20. Marezio, M.; McWhan, D.; Dernier, P.; Remeika, J. Structural aspects of the metal-insulatortransitions in Ti O . Journal of Solid State Chemistry , , 213 – 221.21. Page, Y. L.; Marezio, M. Structural chemistry of magnéli phases Ti n O n − (4 ≤ n ≤ Ti O at 140 K. Journal of Solid State Chemistry , , 13 – 21.22. Lakkis, S.; Schlenker, C.; Chakraverty, B.; Buder, R.; Marezio, M. Metal-insulator transitionsin Ti O single crystals: Crystal characterization, specific heat, and electron paramagneticresonance. Physical Review B , , 1429–1440.23. Ceperley, D. M.; Alder, B. J. Ground State of the Electron Gas by a Stochastic Method. Phys.Rev. Lett. , , 566–569.24. Liborio, L.; Harrison, N. Thermodynamics of oxygen defective Magnéli phases in rutile: Afirst-principles study. Physical Review B , , 104104.25. Eyert, V.; Schwingenschlà ˝ugl, U.; Eckern, U. Charge order, orbital order, and electron local-ization in the Magnéli phase Ti O . Chemical Physics Letters , , 151 – 156.26. Liborio, L.; Mallia, G.; Harrison, N. Electronic structure of the Ti O Magnéli phase.
PhysicalReview B , , 245133.27. Weissmann, M.; Weht, R. Electronic and magnetic properties of the different phases of Ti O from density functional theory. Phys. Rev. B , , 144419.28. Zhong, X.; Rungger, I.; Zapol, P.; Heinonen, O. Electronic and magnetic properties ofTi O predicted by self-interaction-corrected density functional theory. Phys. Rev. B , , 115143. 289. Soler, J. M.; Artacho, E.; Gale, J. D.; Garcà a, A.; Junquera, J.; Ordejøsn, P.; Sà ˛anchez-Portal, D. The SIESTA method for ab initio order- N materials simulation. Journal of Physics:Condensed Matter , , 2745.30. Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulombpotential. The Journal of Chemical Physics , , 8207–8215.31. Heyd, J.; Scuseria, G. E. Efficient hybrid density functional calculations in solids: Assessmentof the Heydâ ˘A ¸SScuseriaâ ˘A ¸SErnzerhof screened Coulomb hybrid functional. The Journal ofChemical Physics , , 1187–1192.32. Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.;Fiolhais, C. Atoms, molecules, solids, and surfaces: Applications of the generalized gradientapproximation for exchange and correlation. Phys. Rev. B , , 6671–6687.33. Kasinathan, D.; Kuneš, J.; Koepernik, K.; Diaconu, C. V.; Martin, R. L.; Prodan, I. m. c. D.;Scuseria, G. E.; Spaldin, N.; Petit, L.; Schulthess, T. C. et al. Mott transition of MnO underpressure: A comparison of correlated band theories. Phys. Rev. B , , 195110.34. Anisimov, V. I.; Gunnarsson, O. Density-functional calculation of effective Coulomb interac-tions in metals. Phys. Rev. B , , 7570–7574.35. Anisimov, V. I.; Zaanen, J.; Andersen, O. K. Band theory and Mott insulators: Hubbard U instead of Stoner I . Phys. Rev. B , , 943–954.36. Anisimov, V. I.; Solovyev, I. V.; Korotin, M. A.; Czy˙zyk, M. T.; Sawatzky, G. A. Density-functional theory and NiO photoemission spectra. Phys. Rev. B , , 16929–16934.37. Anisimov, V. I.; Zaanen, J.; Andersen, O. K. Band theory and Mott insulators: Hubbard U instead of Stoner I . Phys. Rev. B , , 943–954.38. Cococcioni, M.; de Gironcoli, S. Linear response approach to the calculation of the effectiveinteraction parameters in the LDA + U method.
Phys. Rev. B , , 035105.299. Perdew, J. P.; Levy, M. Comment on “Significance of the highest occupied Kohn-Sham eigen-value”. Phys. Rev. B , , 16021–16028.40. Filippetti, A.; Spaldin, N. A. Self-interaction-corrected pseudopotential scheme for magneticand strongly-correlated systems. Phys. Rev. B , , 125109.41. Svane, A.; Gunnarsson, O. Transition-metal oxides in the self-interaction corrected density-functional formalism. Phys. Rev. Lett. , , 1148–1151.42. Szotek, Z.; Temmerman, W. M.; Winter, H. Application of the self-interaction correction totransition-metal oxides. Phys. Rev. B , , 4029–4032.43. Ködderitzsch, D.; Hergert, W.; Temmerman, W. M.; Szotek, Z.; Ernst, A.; Winter, H. Ex-change interactions in NiO and at the NiO(100) surface. Phys. Rev. B , , 064434.44. Fischer, G.; Däne, M.; Ernst, A.; Bruno, P.; Lüders, M.; Szotek, Z.; Temmerman, W.; Herg-ert, W. Exchange coupling in transition metal monoxides: Electronic structure calculations. Phys. Rev. B , , 014408.45. Däne, M.; Lüders, M.; Ernst, A.; Ködderitzsch, D.; Temmerman, W. M.; Szotek, Z.; Herg-ert, W. Self-interaction correction in multiple scattering theory: application to transition metaloxides. Journal of Physics: Condensed Matter , , 045604.46. Hughes, I. D.; DÃd’ne, M.; Ernst, A.; Hergert, W.; LÃijders, M.; Staunton, J. B.; Szotek, Z.;Temmerman, W. M. Onset of magnetic order in strongly-correlated systems from ab ini-tio electronic structure calculations: application to transition metal oxides. New Journal ofPhysics , , 063010.47. Archer, T.; Pemmaraju, C. D.; Sanvito, S.; Franchini, C.; He, J.; Filippetti, A.; Delugas, P.;Puggioni, D.; Fiorentini, V.; Tiwari, R. et al. Exchange interactions and magnetic phases oftransition metal oxides: Benchmarking advanced ab initio methods. Phys. Rev. B , ,115114. 308. Georges, A.; Kotliar, G.; Krauth, W.; Rozenberg, M. J. Dynamical mean-field theory ofstrongly correlated fermion systems and the limit of infinite dimensions. Rev. Mod. Phys. , , 13–125.49. Stanescu, T. D.; Kotliar, G. Fermi arcs and hidden zeros of the Green function in the pseudogapstate. Phys. Rev. B , , 125110.50. Haule, K. Quantum Monte Carlo impurity solver for cluster dynamical mean-field theory andelectronic structure calculations with adjustable cluster base. Phys. Rev. B , , 155113.51. White, S. R. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett. , , 2863–2866.52. Foulkes, W. M. C.; Mitas, L.; Needs, R. J.; Rajagopal, G. Quantum Monte Carlo simulationsof solids. Rev. Mod. Phys. , , 33–83.53. Morales, M. A.; McMinis, J.; Clark, B. K.; Kim, J.; Scuseria, G. E. Multideterminant WaveFunctions in Quantum Monte Carlo. Journal of Chemical Theory and Computation , ,2181–2188.54. Benali, A.; Shulenburger, L.; Romero, N. A.; Kim, J.; von Lilienfeld, O. A. Application of Dif-fusion Monte Carlo to Materials Dominated by van der Waals Interactions. Journal of Chemi-cal Theory and Computation , , 3417–3422.55. Dubecká, M.; JureÄ ka, P.; Derian, R.; Hobza, P.; Otyepka, M.; Mitas, L. Quantum MonteCarlo Methods Describe Noncovalent Interactions with Subchemical Accuracy. Journal ofChemical Theory and Computation , , 4287–4292, PMID: 26589147.56. Ganesh, P.; Kim, J.; Park, C.; Yoon, M.; Reboredo, F. A.; Kent, P. R. C. Binding and Diffusionof Lithium in Graphite: Quantum Monte Carlo Benchmarks and Validation of van der WaalsDensity Functional Methods. Journal of Chemical Theory and Computation , , 5318–5323, PMID: 26583215. 317. Sorella, S.; Casula, M.; Rocca, D. Weak binding between two aromatic rings: Feeling thevan der Waals attraction by quantum Monte Carlo methods. The Journal of Chemical Physics , .58. Send, R.; Valsson, O.; Filippi, C. Electronic Excitations of Simple Cyanine Dyes: ReconcilingDensity Functional and Wave Function Methods. Journal of Chemical Theory and Computa-tion , , 444–455, PMID: 26596164.59. Coccia, E.; Chernomor, O.; Barborini, M.; Sorella, S.; Guidoni, L. Molecular Electrical Prop-erties from Quantum Monte Carlo Calculations: Application to Ethyne. Journal of ChemicalTheory and Computation , , 1952–1962, PMID: 26593830.60. Shulenburger, L.; Mattsson, T. R. Quantum Monte Carlo applied to solids. Phys. Rev. B , , 245117.61. Hood, R. Q.; Kent, P. R. C.; Reboredo, F. A. Diffusion quantum Monte Carlo study of theequation of state and point defects in aluminum. Phys. Rev. B , , 134109.62. Santana, J. A.; Krogel, J. T.; Kim, J.; Kent, P. R. C.; Reboredo, F. A. Structural stabilityand defect energetics of ZnO from diffusion quantum Monte Carlo. The Journal of ChemicalPhysics , , –.63. Wagner, L. K. Ground state of doped cuprates from first-principles quantum Monte Carlocalculations. Phys. Rev. B , , 161116.64. Foyevtsova, K.; Krogel, J. T.; Kim, J.; Kent, P. R. C.; Dagotto, E.; Reboredo, F. A. Ab initio
Quantum Monte Carlo Calculations of Spin Superexchange in Cuprates: The BenchmarkingCase of Ca CuO . Phys. Rev. X , , 031003.65. Devaux, N.; Casula, M.; Decremps, F.; Sorella, S. Electronic origin of the volume collapse incerium. Phys. Rev. B , , 081101. 326. Mostaani, E.; Drummond, N. D.; Fal’ko, V. I. Quantum Monte Carlo Calculation of the Bind-ing Energy of Bilayer Graphene. Phys. Rev. Lett. , , 115501.67. Shulenburger, L.; Baczewski, A.; Zhu, Z.; Guan, J.; Tomà ˇC ˛anek, D. The Nature of theInterlayer Interaction in Bulk and Few-Layer Phosphorus. Nano Letters , , 8170–8175,PMID: 26523860.68. Anderson, J. B. Quantum chemistry by random walk: Higher accuracy. The Journal of Chem-ical Physics , , 3897–3899.69. Umrigar, C. J.; Nightingale, M. P.; Runge, K. J. A diffusion Monte Carlo algorithm with verysmall time-step errors. J. Chem. Phys. , , 2865–2890.70. Hammond, B.; Lester, W.; Reynolds, P. Monte Carlo Methods in Ab Initio Quantum Chem-istry ; World Scientiï ˇn ˛Ac, London, 1994.71. Kim, J.; Esler, K. P.; McMinis, J.; Morales, M. A.; Clark, B. K.; Shulenburger, L.; Ceper-ley, . D. M. Hybrid algorithms in quantum Monte Carlo.
Journal of Physics: ConferenceSeries , , 012008.72. Kim, J.; Esler, K.; McMinis, J.; Ceperley, D. M. Quantum Monte Carlo algorithms: makingmost of large-scale multi/many-core clusters. Chattanooga, TN, 2010.73. Esler, K.; Kim, J.; Shulenburger, L.; Ceperley, D. Fully accelerating quantum Monte Carlosimulations of real materials on GPU clusters. Computing in Science and Engineering , , 40–51.74. Giannozzi et al, P. J. of Physics , , 395502.75. Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations formany-electron systems. Phys. Rev. B , , 5048–5079.76. Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B , , 5188–5192. 337. Umrigar, C. J.; Toulouse, J.; Filippi, C.; Sorella, S.; Hennig, R. G. Alleviation of the Fermion-Sign Problem by Optimization of Many-Body Wave Functions. Phys. Rev. Lett. , ,110201.78. Mitáš, L.; Shirley, E. L.; Ceperley, D. M. Nonlocal pseudopotentials and diffusion MonteCarlo. jcp , , 3467–3475.79. Hamann, D. R.; Schlüter, M.; Chiang, C. Norm-Conserving Pseudopotentials. Phys. Rev. Lett. , , 1494–1497.80. Walter, E. J. Opium Pseuopotential Generation Project. http://opium.sourceforge.net.81. Krogel, J. T.; Santana, J. A.; Reboredo, F. A. Pseudopotentials for quantum Monte Carlostudies of transition metal oxides. Phys. Rev. B , , 075143.82. Rappe, A. M.; Rabe, K. M.; Kaxiras, E.; Joannopoulos, J. D. Optimized pseudopotentials. Phys. Rev. B , , 1227–1230.83. Williamson, A. J.; Hood, R. Q.; Grossman, J. C. Linear-Scaling Quantum Monte Carlo Cal-culations. Phys. Rev. Lett. , , 246406.84. Ewald, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annalen derPhysik , , 253–287.85. Fraser, L. M.; Foulkes, W. M. C.; Rajagopal, G.; Needs, R. J.; Kenny, S. D.; Williamson, A. J.Finite-size effects and Coulomb interactions in quantum Monte Carlo calculations for homo-geneous systems with periodic boundary conditions. Phys. Rev. B , , 1814–1832.86. Chiesa, S.; Ceperley, D. M.; Martin, R. M.; Holzmann, M. Finite-Size Error in Many-BodySimulations with Long-Range Interactions. Phys. Rev. Lett. , , 076404.87. Lin, C.; Zong, F. H.; Ceperley, D. M. Twist-averaged boundary conditions in continuum quan-tum Monte Carlo algorithms. Phys. Rev. E , , 016702.348. Moussa, J. E.; Schultz, P. A.; Chelikowsky, J. R. Analysis of the Heyd-Scuseria-Ernzerhofdensity functional parameter space. The Journal of Chemical Physics , .89. Liechtenstein, A. I.; Anisimov, V. I.; Zaanen, J. Density-functional theory and strong interac-tions: Orbital ordering in Mott-Hubbard insulators. Phys. Rev. B , , R5467–R5470.90. Shick, A. B.; Liechtenstein, A. I.; Pickett, W. E. Implementation of the LDA+U method us-ing the full-potential linearized augmented plane-wave basis. Phys. Rev. B , , 10763–10769.91. Lee, K.-W.; Kuneš, J.; Novak, P.; Pickett, W. E. Disproportionation, Metal-Insulator Transi-tion, and Critical Interaction Strength in Na / CoO . Phys. Rev. Lett. , , 026403.92. Pemmaraju, C. D.; Archer, T.; Sánchez-Portal, D.; Sanvito, S. Atomic-orbital-based approx-imate self-interaction correction scheme for molecules and solids. Phys. Rev. B , ,045101.93. Akande, A.; Sanvito, S. Exchange parameters from approximate self-interaction correctionscheme. The Journal of Chemical Physics ,127