Free energy from stationary implementation of the DFT+DMFT functional

1.1 MB·English

by  Kristjan Haule

Free energy from stationary implementation of the DFT+DMFT functional
Kristjan Haule and Turan Birol
Department of Physics and Astronomy, Rutgers University, Piscataway, USA
(Dated: January 29, 2015)
The stationary functional of the all-electron density functional plus dynamical mean field theory
(DFT+DMFT) formalism to perform free energy calculations and structural relaxations is imple-
mented for the first time. Here, the first order error in the density leads to a much smaller, second
order error in the free energy. The method is applied to several well known correlated materials;
metallic SrVO , Mott insulating FeO, and elemental Cerium, to show that it predicts the lattice
3
constantswithveryhighaccuracy. InCerium,weshowthatourmethodpredictstheiso-structural
transition between the α and γ phases, and resolve the long standing controversy in the driving
5 mechanism of this transiton.
1
0 PACSnumbers: 71.27.+a,71.30.+h
2

n Prediction of the crystal structures of solids by large which delivers stationary free energies at finite temper-
a
scale quantum mechanical simulations is one of the fun- atures. This stationarity is crucial for practical imple-
J
damentalproblemsofcondensedmatterphysics, andoc- mentation and precision of computed energies, since the
7
2 cupiesacentralplaceinmaterialsdesign. Theworkhorse firstordererrorinthedensityρ(ortheGreen’sfunction)
ofthefieldistheDensityFunctionalTheory(DFT)[1]at leads only to the much smaller second order error in the

] thelevelofLocalDensityApproximation(LDA)orGen- free energy, since the first order variation vanishes, i.e.,
l
e eralized Gradient Approximations (GGAs), which pre- δF/δρ = 0. This property is also crucial in calculating
-
r dict lattice constants of weakly correlated materials typ- the forces, as stationarity of the functional ensures that
t ically within 1% relative error [2]. only Hellmann-Feynman forces need to be computed for
s
∼
. structural relaxation . [58]
t These errors of DFT in LDA/GGA implementations
a The DFT+DMFT total energy is given by [9] :
are an order of magnitude larger in the so called cor-
m
related materials: For example, the lattice constant of
- 1
d δ-Pu is underestimated by 11% [3] or non-magnetic FeO E =Tr(H0G)+ Tr(ΣG)+EH[ρ]+Exc[ρ]
2
n by7%[4]. WhileGGAsandhybridfunctionalscansome-
o times improve upon conventional LDA, these function- −ΦDC[nloc]+Enuc−nuc (1)
c
alsmanytimesdegradetheagreementbetweenpredicted
[ where H = 2 + δ(r r(cid:48))V (r), G is the elec-
1 aconndsteaxnptesr,iminepnatartlliycudlaerteirnmminaetderibaulslkcomnotadiunliinagnhdealavtytiecle- tron Gre0en’s f−un∇ction, EH−[ρ] anedxtExc[ρ] are Hartree
v ements. [2] and DFT exchange-correlation functional, Vext is the
6 electron-nuclearpotential,E istheinteractionen-
nuc−nuc
3 To account for the correlation effects, more sophisti- ergyofnuclei,ΣistheDMFTself-energy,andΦDC[n ]
loc
9 cated many body methods have been developed. Among
is the double-counting (DC) functional. [4] Here the
6
them, one of the most successful algorithms is the dy- Migdal-Galitskiiformula(MGF)isusedE = 1Tr(ΣG)
0 pot 2
namical mean-field theory (DMFT) [5]. It replaces the
. to compute the DMFT part of the potential energy.
1 problem of describing correlation effects in a periodic
Gordon Baym showed [31] that for certain class of ap-
0
lattice by a strongly interacting impurity coupled to a
5 proximations, which are derivable from a functional ex-
self-consistent bath [6]. To become material specific,
1 pressedintermsofclosed-loopFeynmandiagrams,MGF
: DMFT was soon developed into an electronic structure can be used instead of more complicated expression for
v
tool (LDA+DMFT) [7, 8], which achieved great success
i evaluating the Luttinger-Ward Functional [32, 33]. He
X in numerous correlated materials (for a review see [9]).
calledsuchapproximationsconserving. WhiletheDMFT
r The LDA+DMFT method has mainly been used for the is a conserving approximation in Baym’s sense, LDA or
a
calculation of spectroscopic quantities, and only a few GGA are not, as the Galitskii-Migdal formula 1Tr(V ρ)
2 xc
dozens[10–30]ofstudiesmanagedtocomputeenergetics
hastobereplacedbytheexchange-correlationfunctional
of correlated solids, and only a handful of them used ex- Exc[ρ]. As a result, the combination of DFT+DMFT in
actsolversandchargeself-consistency[18,19,24,25,28,
itscharge-selfconsistentversionisnotconservingeither,
29]. This is not only because of the very high computa-
and consequently MGF can give different total energy
tionalcost,butalsobecausepreviousimplementationsof
thantheLuttinger-Wardfunctional. Onlytheevaluation
LDA+DMFTwerenotstationary,andhenceitwashard
of the latter is guaranteed to give stationary free ener-
to achieve precision of free energies needed for structure
gies. We will give numerical evidence that evaluation of
optimization and study of phase transitions in solids.
MGF in Eq. 1 gives different results than evaluation of
Here we implemented the LDA+DMFT functional, the Luttinger-Ward functional, which strongly suggests
2
that Eq. 1 gives non-stationary total energies. where Z can be directly computed from atomic en-
atom
The Luttinger-Ward functional of DFT+DMFT has ergies, and P is the probability for no kinks on the im-
0
been well known for several years [9], but it has never purity. Finally, Z =exp( F /T), giving directly
imp
−
been successfully implemented to compute the free en-
ergy of a solids. It has the following form Fimp = T(log(Zatom) log(P0)). (5)
− −
Γ[G]=TrlogG Tr((G−1 G−1)G)+EH[ρ] When the temperature is low, P0 becomes exponen-
− 0 − tially small, and we can no longer determine Z to high
+Exc[ρ]+ΦDMFT[PˆG] ΦDC[Pˆρ]+E , (2)
− nuc−nuc enough precision in this way. However, we can compute
very precisely the internal energy of the impurity at ar-
where G−1(rr(cid:48);iω) = [iω + µ + 2 V (r)]δ(r
0 ∇ − ext − bitrary temperature. The internal energy of QIM E
r(cid:48)), ΦDMFT[PˆG] is the DMFT functional, which is imp
is given by
the sum of all local skeleton Feynman diagrams.
The projected Green’s function PˆG Glocal = d∆
(cid:80) ≡ E =Tr((∆+ε ω )G )+E , (6)
LL(cid:48)|φL(cid:105)(cid:104)φL|G|φL(cid:48)(cid:105)(cid:104)φL(cid:48)| and the projected density imp imp− ndωn imp imp−pot
Pˆρ ρ are computed with projection to a set of lo-
local
caliz≡ed functions φ centered on the ”correlated” atom. which follows directly from the thermodynamic av-
The projection de|fin(cid:105)es the local Green’s function G , erage of QIM Hamiltonian. Here the hybridization
local
the essential variable of the DMFT. ∆ and impurity levels εimp are determined from the
The variation of functional Γ[G] with respect to G local green’s function by the standard DMFT self-
(δΓ[G]/δG) gives, consistency condition G−lo1cal = iωn − εimp − Σ − ∆,
and E = 1Tr(ΣG ). These quantities can be
imp−pot 2 imp
G−1 G−1+(V +V )δ(r r(cid:48))δ(τ τ(cid:48)) computed very precisely by CTQMC using the follow-
− 0 H xc − −
ing tricks: i) Tr(∆G ) is computed from the average
δΦDMFT[G ] imp
+Pˆ local perturbation order k of CTQMC, and takes the form
δG (cid:104) (cid:105)
local Tr(∆G ) = k /T, where T is temperature [34]; ii)
imp
δΦDC[ρ ] (cid:104) (cid:105)
Pˆ local δ(r r(cid:48))δ(τ τ(cid:48))=0, (3) Eimp−pot is computed from the energies of atomic state
− δρlocal − − of QIM Ematom and their probabilities Pm by Eimp−pot =
(cid:80) P Eatom Tr(ε n ) [34], which delivers much
which vanishes, since it is equal to the Dyson equation m m m − imp imp
more precise interaction energy than obtained by MGF;
thatdeterminesself-consistentG,hencethefunctionalis
ii) We spline ∆(ω ) in Matsubara points and determine
n
stationary.
its derivative d∆/dω , and then carry out Matsubara
n
The value of the functional Γ at the self-consistently
sum by subtracting out the leading high-frequency tails
determined G delivers the free energy of the system [31].
by formula A/((iω ε )(iω ε )), which hasan analytic
We evaluate it by inserting G−1 G−1 from Eq. 3 into − 1 − 2
0 − sum of A(f(ε1) f(ε2))/(ε1 ε2). Because probabilities
Eq. 2 to obtain − −
P and perturbation order k are known to very high
m
(cid:104) (cid:105)
precision in CTQMC, the impurity internal energy can
F =E Tr((V +V )ρ)+EH[ρ]+Exc[ρ]
nuc−nuc− H xc easilybecomputedwithprecisionofafractionofameV.
+TrlogG−TrlogGloc+Fimp TocomputepreciseimpurityfreeenergyFimp atlower
+Tr(V ρ ) ΦDC[ρ ]+µN, (4) temperature, we first converge DFT+DMFT equations
dc loc loc
−
to high accuracy at low temperature. Using converged
where we denoted Vdc ≡ δΦDC[ρlocal]/δρlocal and Fimp impurityhybridization∆(iωn),weraisethetemperature
is the free energy of the impurity problem, i.e., Fimp = of the impurity (keeping ∆ fixed) to T>, so that P0 is
TrlogGloc−Tr(ΣGloc)+ΦDMFT[Gloc]. [4] Here we also of the order of 10−5 or higher, and obtain reliable Fimp
use the fact the solution of the auxiliary impurity prob- using Eq. 5, and entropy S at this higher tempera-
imp
lem delivers the exact local Green’s function, i.e., Σ = tureS =(E (T ) F (T ))/T . Next, weevalute
> imp > imp > >
δΦDMFT[Glocal]/δGlocal, and we added µN because we impurityinternalener−gyforseveralinversetemperatures
work at constant electron number. β = 1/T, and than we use standard thermodynamic re-
Thecrucialpointisthatthecontinuoustimequantum lations to obtain entropy at lower temperature T by
MonteCarlomethod(CTQMC)[34,35]solvesthequan-
tum impurity model (QIM) numerically exactly, hence, E (T ) E (T) (cid:90) 1/T
imp > imp
S(T)=S + E (β)dβ(7)
we can compute very precisely the impurity internal en- >− T T − imp
> 1/T>
ergyaswellasthefreeenergyF ofthismodel. Athigh
imp
enough temperature, F can be directly read-off from where β = 1/T. This formula is obtained integrating
imp
(cid:82)
the probability for the perturbation order k, which we by parts the standard formula S = c /TdT and c =
v v
call P , and is computed by P =Z /Z (where Z is the dE/dT. WehenceobtainS andF =E TS
k k k k imp imp imp imp
−
partition function with k kinks), hence Z = Z /P , at lower T which can be inserted into Eq. 4. The rest
atom 0
−
3
of the terms in Eq. 4 are relatively straightforward to
evaluate, however, for a high precision implementation
one needs to combine the terms that largely cancel and
evaluate them together [4].
Previous implementations of free energy within
LDA+DMFT [25, 27, 36] were based on i) evaluating
the total energy Eq. 1 at range of temperatures, and in-
tegratingresultingspecificheat[36], andii)thecoupling
constant integration [25, 27], where total energy of the
solid is needed for a range of coulomb repulsion’s U and
is than integrated over U. In both approaches, the self-
consistentLDA+DMFTsolutionisneededformanyval-
ues of the parameters (either U or T) to evaluate F. In
our method, a single LDA+DMFT calculation for solid
is needed, which makes the method much more efficient.
Furthermore, current implementation of the free energy
is stationary, hence higher precision of F is achieved.
To test the implementation of the LDA+DMFT func-
tional, we computed the volume dependence of the free
energy for three well studied correlated materials: a
metallic early transition metal oxide with perovskite
structure SrVO , a Mott insulating transition metal ox-
3
FIG. 1: (Color online): a) E(V) and F(V) for SrVO at
ide FeO in its rock salt structure, and the lanthanide 3
T = 230K from Eq. 1 and 4, respectively. Entropy term
elementalmetal,Cerium,initsfacecenteredcubicstruc-
TSimp(V)isverysmall. (b)theoreticalandexperimental[41]
ture, which undergoes a first order iso-structural transi-
p(V). Good agreement between theoretical −dF/dV and ex-
tion. periment is found. (c) Impurity Entropy Eq. 7 for represen-
We used the implementation of LDA+DMFT of tative volumes. To obtain Simp, temperature is varied in the
Ref.37, which is based on the Wien2K package [38], and impurityproblemonly,andnotintheLDA+DMFTproblem
of the solid.
LDA in combination with nominal double-counting [37,
39, 40]. More technical details are given in the supple-
mentary information.
notice a temperature scale at which t2g shell is degen-
In Fig. 1(a) we show the energy E(V), and F(V) erate (log(6)) nor the scale of the lowest order Kramers
for SrVO3 at T = 230K, computed with Eq. 1, and doublet (log(2)), but we notice the Fermi liquid scale in
Eq. 4, respectively. The minima of E(V) and F(E) the steep downturn of S(T) at T 350K.
are achieved at 55.71˚A3 and 55.51˚A3. The experimen- Fig.2(a)showsE(V)andF(V)≈forparamagneticMott
tally determined volume is V = 56.53˚A3 [42]. The insulating FeO at 300K, above its antiferromagnetic or-
exp
LDA+DMFT hence slightly underestimates the equilib- dering temperature. The equilibrium volume of E and
rium volume (1.8%), which gives 0.6% error in lattice F is 20.28˚A3 and 20.24˚A3, while the experimental vol-
constant. This is well within the standard error of best ume is 20.342˚A3. The lattice constant is thus under-
DFT functionals for weakly correlated materials. estimated for only 0.10% and 0.16%, respectively. In
ThemetallicnatureofSrVO3,withmoderatemassen- comparison, all standard DFT functionals severally un-
hancementsm∗/mband 2[4],leadstoverysmallDMFT derestimate FeO lattice constant, for example PBE-sol,
≈
correctionsincrystalstructure[4]. Notethatenergymin- PBE, and LDA for 5.2%, 5.0%, and 7.7%, respectively.
imizationleadstoslightlylargervolumethanfreeenergy In Fig. 2(b) we show P(V) diagram and its excellent
minimization, contrary to expectations. This is because agreement with experiment. Fig. 2(c) shows impurity
energyiscomputedfromnon-stationaryEq.1,whilefree- entropyS (T)forafewvolumes. Incontrasttometal-
imp
energy is obtained from the stationary expression Eq. 4. lic SrVO , here we clearly see an extended plateau of
3
The latter is hence more trustworthy, and should be S (T)=log(6) k around1000K,whichsignalscom-
imp B
∗
considered best LDA+DMFT result. This is also clear plete degeneracy of the t shell, and its slight decrease
2g
from pressure versus volume diagram in Fig. 1(b), where at 300K in proximity to the AFM state.
dF/dV agreesmorefavorablywiththeexperimentthan The iso-structural transitions of Cerium attracted a
−
dE/dV obtained by MGF. lot of experimental and theoretical effort, but its theo-
−
In Fig. 1(c), we show the impurity entropy obtained retical understanding is still controversial. On the ba-
by Eq. 7 for two representative volumes. In this itin- sis of LDA+DMFT calculation McMahan et.al [11] pro-
erant system with very large hybridization, we do not posed that the total energy exhibits a double-minimum
4
shape,concomitantwiththeappearanceofthequasipar-
ticlepeakattemperatureashighas1500K,signalingthe
first order transition. Using different implementation of
the same method, Amadon et.al [27, 46] proposed that
the transition is entropy driven, and that the total en-
ergy is featureless with the minimum corresponding to
low volume α-phase. Only the addition of the entropy
termmovestheminimumtothelargervolumeofγ-phase.
In this picture the transition at low temperatures, where
the entropy becomes small and cannot drive the tran-
sition, is intrinsically absent. Yet another proposal was
recentlyputforwardonthebasisofLDA+Gutzwillercal-
culations [47, 48], in which the transition is present even
atzerotemperature,butthetransitionoccursatnegative
pressure. Thetransitionisthusdetectableevenintheto-
talenergy,intheabsenceofentropy,andbecomessecond
order at T = 0. In the same method, the finite temper-
ature transition is first order, and the double-minimum
shape of free energy becomes most pronounced at very
high temperature (1500K) [48].
Our LDA+DMFT results for Ce are plotted in Fig. 3.
The total energy curve at 400K clearly shows a region
of very flat shape in the region between the α-γ volume.
FIG. 2: (Color online): a) E(V) and F(V) for FeO from Indeed the derivative of the energy dE/dV displayed
−
Eq. 1 and 4, respectively. Entropy term TSimp(V) is large in Fig. 3(c) shows a clear region of zero slope around
but almost constant. (b) theoretical and experimental p(V). 1GPa. ThisisconsistentwithresultsofLanataetal.[47]
FilledandemptycirclesarefromRefs.43and44,respectively.
finding very similar zero slope of dE/dV at zero tem-
(c) Impurity entropy Eq. 7 for representative volumes. The −
perature, but is inconsistent with Ref. 27, which finds
degeneracy of the t2g shell above 1000K is apparent.
no feature in total energy. It is also inconsistent with
McMahan et.al [11] showing clear double-peak in total
energy. On the other hand, the addition of entropy sub-
stantiallyincreasetheregionofsoftvolume,assuggested
by Amadon et.al [46]. Indeed the change of the entropy
between the two phases is of the order of 0.9k , which
B
is consistent with experimental estimations of 30meV at
400K [49]. The physical mechanism behind this large
entropy change and unusual volume dependence of en-
ergy is in very fast variation of coherence temperature,
as suggested in Refs. [11, 46], and conjectured in Kondo
volume collapse theory [50]. The phase transition in our
calculation occurs around 1.6GPa, which is not far from
experimentally determined critical pressure of 1.25GPa
at T =400K. The free energy barrier in our calculation
is however extremely small, and no clear double peak
of F(V) or negative slope of dF/dV can be detected
−
withinour1meVprecisionofenergies. Thisissimilarto
results of Ref. 48 at 400K, but different from Ref. 11.
While the start of the transition region in α-phase is
in good agreement with experiment, the γ-phase vol-
ume is underestimated in our calculation. We believe
that the addition of phonon entropy is needed to further
FIG. 3: (Color online): a) E(V) and F(V) for elemental increase the transition region, and establish larger free
CeriumfromEq.1and 4,respectively. Dataarepresentedfor energy barrier between the two phases. Experimentally,
T=400 and 900K. (b) Entropy Simp(V) is large and changes above 460K the α γ phase transition ends with the fi-
dramaticallyaccrosthetransiton. (c)theoreticalandexperi- nite temperature c−ritical point. Our calculation at high
mental [45] p(V) diagram.
temperature900Kshowsthatthesignatureofthephase
5
transition in F(V) and E(V) disappears, which is dif- [13] B. Amadon, S. Biermann, A. Georges, and F. Aryaseti-
ferent than predicted by Gutzwiller method [48], where awan,Phys.Rev.Lett.96,066402(2006),URLhttp://
the largest free energy barrier is found at these elevated link.aps.org/doi/10.1103/PhysRevLett.96.066402.
[14] S. Y. Savrasov, K. Haule, and G. Kotliar, Phys. Rev.
temperatures, but qualitatively consistent with Ref. 11 .
Lett. 96, 036404 (2006), URL http://link.aps.org/
In summary, we successfully implemented the station-
doi/10.1103/PhysRevLett.96.036404.
ary formula for the free energy of DFT+DMFT method. [15] L. V. Pourovskii, B. Amadon, S. Biermann, and
OntheexampleofSrVO3,FeOandCemetalwedemon- A.Georges,Phys.Rev.B76,235101(2007),URLhttp:
strated that the method successfully predicts lattice vol- //link.aps.org/doi/10.1103/PhysRevB.76.235101.
umesincorrelatedsolids,whicharedifficultforstandard [16] I. Leonov, N. Binggeli, D. Korotin, V. I. Anisi-
mov, N. Stoji´c, and D. Vollhardt, Phys. Rev. Lett.
DFT functionals. We also resolved controversy in the
101, 096405 (2008), URL http://link.aps.org/doi/
mechanism of the α-γ transition in Cerium.
10.1103/PhysRevLett.101.096405.
This work was supported Simons foundation under
[17] I. Di Marco, J. Min´ar, S. Chadov, M. I. Katsnel-
project ”Many Electron Problem”, and by NSF-DMR son, H. Ebert, and A. I. Lichtenstein, Phys. Rev. B
1405303. T.B. was supported by the Rutgers Center for 79,115111(2009),URLhttp://link.aps.org/doi/10.
Materials Theory. This research used resources of the 1103/PhysRevB.79.115111.
Oak Ridge Leadership Computing Facility at the Oak [18] M.Aichhorn, L.Pourovskii, andA.Georges, Phys.Rev.
B 84, 054529 (2011), URL http://link.aps.org/doi/
Ridge National Laboratory, which is supported by the
10.1103/PhysRevB.84.054529.
Office of Science of the US Department of Energy under
[19] G. Lee, H. S. Ji, Y. Kim, C. Kim, K. Haule, G. Kotliar,
Contract No. DE-AC05-00OR22725.
B. Lee, S. Khim, K. H. Kim, K. S. Kim, et al., Phys.
Rev. Lett. 109, 177001 (2012), URL http://link.aps.
org/doi/10.1103/PhysRevLett.109.177001.
[20] B. Amadon, Journal of Physics: Condensed Mat-
ter 24, 075604 (2012), URL http://stacks.iop.org/
[1] P. Hohenberg and W. Kohn, Phys. Rev. 136, 0953-8984/24/i=7/a=075604.
B864(1964),URLhttp://link.aps.org/doi/10.1103/ [21] I. Leonov, A. I. Poteryaev, V. I. Anisimov, and D. Voll-
PhysRev.136.B864. hardt, Phys. Rev. B 85, 020401 (2012), URL http:
[2] A. E. Mattsson, R. Armiento, J. Paier, G. Kresse, //link.aps.org/doi/10.1103/PhysRevB.85.020401.
J. M. Wills, and T. R. Mattsson, The Journal [22] M.S.Litsarev,I.DiMarco,P.Thunstro¨m,andO.Eriks-
of Chemical Physics 128, 084714 (2008), URL son,Phys.Rev.B86,115116(2012),URLhttp://link.
http://scitation.aip.org/content/aip/journal/ aps.org/doi/10.1103/PhysRevB.86.115116.
jcp/128/8/10.1063/1.2835596. [23] O. Grns, I. D. Marco, P. Thunstrm, L. Nordstrm,
[3] S. Y. Savrasov, G. Kotliar, and E. Abrahams, Nature O. Eriksson, T. Bjrkman, and J. Wills, Computa-
410, 793 (2001), URL http://dx.doi.org/10.1038/ tional Materials Science 55, 295 (2012), ISSN 0927-
35071035. 0256, URL http://www.sciencedirect.com/science/
[4] For more information see the supplementary informa- article/pii/S092702561100646X.
tion. [24] D. Grieger, C. Piefke, O. E. Peil, and F. Lechermann,
[5] A. Georges and G. Kotliar, Phys. Rev. B 45, Phys.Rev.B86,155121(2012),URLhttp://link.aps.
6479 (1992), URL http://link.aps.org/doi/10.1103/ org/doi/10.1103/PhysRevB.86.155121.
PhysRevB.45.6479. [25] L. V. Pourovskii, T. Miyake, S. I. Simak, A. V. Ruban,
[6] A.Georges,G.Kotliar,W.Krauth,andM.J.Rozenberg, L. Dubrovinsky, and I. A. Abrikosov, Phys. Rev. B
Rev.Mod.Phys.68,13(1996),URLhttp://link.aps. 87,115130(2013),URLhttp://link.aps.org/doi/10.
org/doi/10.1103/RevModPhys.68.13. 1103/PhysRevB.87.115130.
[7] V. I. Anisimov, A. I. Poteryaev, M. A. Korotin, A. O. [26] I. Leonov, V. I. Anisimov, and D. Vollhardt, Phys. Rev.
Anokhin,andG.Kotliar,JournalofPhysics: Condensed Lett. 112, 146401 (2014), URL http://link.aps.org/
Matter 9, 7359 (1997), URL http://stacks.iop.org/ doi/10.1103/PhysRevLett.112.146401.
0953-8984/9/i=35/a=010. [27] J. Bieder and B. Amadon, Phys. Rev. B 89,
[8] A.I.LichtensteinandM.I.Katsnelson,PhysicalReview 195132 (2014), URL http://link.aps.org/doi/10.
B 57, 6884 (1998), URL http://dx.doi.org/10.1103/ 1103/PhysRevB.89.195132.
PhysRevB.57.6884. [28] S. Mandal, R. E. Cohen, and K. Haule, Phys. Rev. B
[9] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, 89,220502(2014),URLhttp://link.aps.org/doi/10.
O.Parcollet,andC.A.Marianetti,Rev.Mod.Phys.78, 1103/PhysRevB.89.220502.
865 (2006), URL http://link.aps.org/doi/10.1103/ [29] H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev.
RevModPhys.78.865. B 89, 245133 (2014), URL http://link.aps.org/doi/
[10] S. Y. Savrasov, G. Kotliar, and E. Abrahams, Nature 10.1103/PhysRevB.89.245133.
410, 793 (2001). [30] H. Park, A. J. Millis, and C. A. Marianetti, Phys. Rev.
[11] K. Held, A. K. McMahan, and R. T. Scalettar, Phys. B 90, 235103 (2014), URL http://link.aps.org/doi/
Rev. Lett. 87, 276404 (2001), URL http://link.aps. 10.1103/PhysRevB.90.235103.
org/doi/10.1103/PhysRevLett.87.276404. [31] G. Baym, Phys. Rev. 127, 1391 (1962), URL http://
[12] A. K. McMahan, K. Held, and R. T. Scalettar, Phys. link.aps.org/doi/10.1103/PhysRev.127.1391.
Rev. B 67, 075108 (2003), URL http://link.aps.org/ [32] J. M. Luttinger and J. C. Ward, Phys. Rev. 118,
doi/10.1103/PhysRevB.67.075108. 1417 (1960), URL http://link.aps.org/doi/10.1103/
6
PhysRev.118.1417. 026403 (2012), URL http://link.aps.org/doi/10.
[33] G. Baym and L. P. Kadanoff, Phys. Rev. 124, 1103/PhysRevLett.108.026403.
287 (1961), URL http://link.aps.org/doi/10.1103/ [52] A. McMahan, C. Huscroft, R. Scalettar, and E. Pol-
PhysRev.124.287. lock,JournalofComputer-AidedMaterialsDesign5,131
[34] K. Haule, Phys. Rev. B 75, 155113 (2007), URL http: (1998), ISSN 0928-1045, URL http://dx.doi.org/10.
//link.aps.org/doi/10.1103/PhysRevB.75.155113. 1023/A%3A1008698422183.
[35] P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and [53] M.Weinert,E.Wimmer,andA.J.Freeman,Phys.Rev.
A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006), URL B26,4571(1982),URLhttp://link.aps.org/doi/10.
http://link.aps.org/doi/10.1103/PhysRevLett.97. 1103/PhysRevB.26.4571.
076405. [54] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov,
[36] K. Held, A. K. McMahan, and R. T. Scalettar, Phys. M. Troyer, and P. Werner, Rev. Mod. Phys. 83,
Rev. Lett. 87, 276404 (2001), URL http://link.aps. 349 (2011), URL http://link.aps.org/doi/10.1103/
org/doi/10.1103/PhysRevLett.87.276404. RevModPhys.83.349.
[37] K. Haule, C.-H. Yee, and K. Kim, Phys. Rev. B [55] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.
81,195107(2010),URLhttp://link.aps.org/doi/10. Lett. 77, 3865 (1996).
1103/PhysRevB.81.195107. [56] J.P.Perdew,A.Ruzsinszky,G.I.Csonka,O.A.Vydrov,
[38] P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, G.E.Scuseria,L.A.Constantin,X.Zhou,andK.Burke,
and J. Luitz, WIEN2K, An Augmented Plane Wave + Phys.Rev.Lett.100,136406(2008),URLhttp://link.
LocalOrbitalsProgramforCalculatingCrystalProperties aps.org/doi/10.1103/PhysRevLett.100.136406.
(Karlheinz Schwarz, Techn. Universit¨at Wien, Austria, [57] M. Takizawa, M. Minohara, H. Kumigashira, D. Toy-
2001). ota, M. Oshima, H. Wadati, T. Yoshida, A. Fuji-
[39] K. Haule, T. Birol, and G. Kotliar, Phys. Rev. B mori, M. Lippmaa, M. Kawasaki, et al., Phys. Rev. B
90,075136(2014),URLhttp://link.aps.org/doi/10. 80,235104(2009),URLhttp://link.aps.org/doi/10.
1103/PhysRevB.90.075136. 1103/PhysRevB.80.235104.
[40] K.Haule,eprintarXiv:1501.03438v1(2015),URLhttp: [58] NotethatPulayforcesappearwhenincompletebasisset
//arxiv.org/abs/1501.03438. is used.
[41] R. Hrubiak, Ph.D. thesis, Florida International Uni-
versity (2012), URL http://digitalcommons.fiu.edu/
cgi/viewcontent.cgi?article=1802&context=etd.
[42] M. Goodenough, J. B.and Longo, SpringerMaterials -
The Landolt-Brnstein Database 4a (1970), URL http:
//dx.doi.org/10.1007/10201420_50.
[43] T. Yagi, T. Suzuki, and S.-I. Akimoto, Journal of
Geophysical Research: Solid Earth 90, 8784 (1985),
ISSN 2156-2202, URL http://dx.doi.org/10.1029/
JB090iB10p08784.
[44] R. L. Clendenen and H. G. Drickamer, The Jour-
nal of Chemical Physics 44, 4223 (1966), URL
http://scitation.aip.org/content/aip/journal/
jcp/44/11/10.1063/1.1726610.
[45] M. J. Lipp, D. Jackson, H. Cynn, C. Aracne, W. J.
Evans, and A. K. McMahan, Phys. Rev. Lett. 101,
165703 (2008), URL http://link.aps.org/doi/10.
1103/PhysRevLett.101.165703.
[46] B. Amadon, S. Biermann, A. Georges, and F. Aryaseti-
awan,Phys.Rev.Lett.96,066402(2006),URLhttp://
link.aps.org/doi/10.1103/PhysRevLett.96.066402.
[47] N. Lanata`, Y.-X. Yao, C.-Z. Wang, K.-M. Ho,
J.Schmalian,K.Haule,andG.Kotliar,Phys.Rev.Lett.
111, 196801 (2013), URL http://link.aps.org/doi/
10.1103/PhysRevLett.111.196801.
[48] N. Lanata`, Y.-X. Yao, C.-Z. Wang, K.-M. Ho, and
G.Kotliar,Phys.Rev.B90,161104(2014),URLhttp:
//link.aps.org/doi/10.1103/PhysRevB.90.161104.
[49] F. Decremps, L. Belhadi, D. L. Farber, K. T. Moore,
F. Occelli, M. Gauthier, A. Polian, D. Antonangeli,
C.M.Aracne-Ruddle,andB.Amadon,Phys.Rev.Lett.
106, 065701 (2011), URL http://link.aps.org/doi/
10.1103/PhysRevLett.106.065701.
[50] J. W. Allen and L. Z. Liu, Phys. Rev. B 46,
5047 (1992), URL http://link.aps.org/doi/10.1103/
PhysRevB.46.5047.
[51] K. Ohta, R. E. Cohen, K. Hirose, K. Haule,
K. Shimizu, and Y. Ohishi, Phys. Rev. Lett. 108,
7
SUPPLEMENTARY INFORMATION: FREE ENERGY FROM STATIONARY IMPLEMENTATION OF
THE DFT+DMFT FUNCTIONAL
TECHNICAL DESCRIPTION OF 3dorbitalsofVanadium. Themuffin-thinradiusofVana-
COMPUTATIONAL DETAILS dium was set to R = 1.83a , and U = 10eV was
mt B
used, which was previously shown to give good spectra
We used the implementation of LDA+DMFT of [39,40]forthislocalizedorbital(seespectrabelow). The
Ref.37, which is based on Wien2K package [38]. The YukawaformofscreeninginteractionthangivesJ 1eV
≈
exchange-correlationfunctionalE ofLDAwasutilized, (see note below). Brillouin zone integrations were done
xc
in combination with nominal double-counting (DC) [37, over 15 15 15 k-point in the whole zone in the self-
× ×
39], which was shown to be closest to the exact form of consistent calculations, and for calculation of the impu-
DC[40]. Wechecked(ontheexampleofCerium)thatthe rity entropy, the hybridization is computed on more pre-
exact-DCgivesverysimilarfreeenergy,asexpectedfora cise 36 36 36 k-points mesh. We mention in pass-
× ×
stationary functional. The convergence of LDA+DMFT ing that impurity entropy is very sensitive to the precise
results is much faster using nominal DC, hence most of frequency dependence of the hybridization, and requires
results in this publication are obtained by this simplifi- very dense momentum mesh.
cation. For FeO, all five 3d orbitals are treated by DMFT and
The impurity model is solved using the hybridization the muffin-thin radius of iron is set to R = 2.11a ,
mt B
expansion version of the numerically exact continuous andtheCoulombrepulsiontopreviouslydeterminedU =
time QMC method [34, 35]. Of the order of 300 LDA 8eV [51], which requires J 1eV in Yukawa form. In
≈
and 30 DMFT iterations were required for precision of Cemetal,allseven4f orbitalsaretreatedbyDMFTand
1meV per formula unit, and between 109 1010 Monte themuffin-thinsphereisR =2.5a ,thek-pointmesh
mt B
−
Carlomoveswereacceptedperimpurityiterationforpre- is 21 21 21, and the Coulomb U = 6eV [27, 46, 52],
× ×
cise enough impurity solution. The resources of Titan leads to J = 0.72eV in Yukawa form. The spin-orbit
supercomputer were used. coupling is included in Cerium, but neglected in SrVO
3
To construct the projector, the atomic-like local or- and FeO.
bitals are used rφ = ul(r)Y (ˆr). The radial part of
(cid:104) | lm(cid:105) r lm
the local orbital u (r) is the solution of the scalar rela-
l
tivistic Dirac equation inside the muffin-tin sphere, lin- DETAILS ON EVALUATION OF LDA+DMFT
earized at the Fermi level. The muffin-tin spheres are FUNCTIONAL
set to touch at the lowest volume. We tested a few other
formsoftheprojectorsdefinedinRef.37. Thestationary HereweexplainhowweevaluatethetotalenergyEq.1
F(V) is quite insensitive to the precise choice of projec- and the free energy Eq.4 of the main text.
tor, however, E(V) changes much more. For total energy Eq.1, we group the terms in the fol-
ForSrVO calculations,wetreateddynamicallyallfive lowing way
3
1
E =Tr(( 2+V +V +V )G) Tr((V +V )ρ)+EH[ρ]+Exc[ρ]+E + Tr(ΣG) ΦDC[ρ ] (8)
ext H xc H xc nuc−nuc loc
−∇ − 2 −
WethensplittheenergyintothreetermsE =E +E + on LDA+DMFT charge. We then evaluate
1 2
E , where the first two parts E , E are computed using
3 1 2 E =Tr(εDFTG) (9)
the Green’s function of the solid, and the third E using 1
3
the impurity Green’s function. and
E = Tr((V +V )ρ)+EH[ρ]+Exc[ρ]+E (.10)
ThefirstfivetermsinEq.8looksimilartothestandard 2 − H xc nu−nuc
DFT energy functional, except that the Green’s func- Both E and E are computed using Green’s function G
1 2
tion G here is the self-consistent LDA+DMFT Green’s anddensityρofthesolidinthesamewayasthestandard
function. We first solve the eigenvalue problem for DFT total energy is implemented [53].
Kohn-Sham states ( 2 + V + V + V )ψ = The last two terms of Eq. 8 can be computed either
ext H xc ik
εDFTψ , where εDFT−a∇re DFT-like energies, computed from the local Green’s function PˆG or from the impu-
ik ik ik
8
rity Green’s function G . Once the self-consistency is thepotentialenergyfromtheimpurityprobabilities,i.e.,
imp
reached, the two are of course equal. We choose to eval-
uate the second term on the impurity Gimp 1 (cid:88)
Tr(Σ G )= P Eatom Tr(ε n )
2 imp imp m m − imp imp
1
E = Tr(Σ G ) ΦDC[ρ ] (11) m
3 imp imp imp
2 −
However,weneveractuallyuseMigdal-Galitskiiformula, The free energy functional Γ[G] (Eq. 2 of the main
becauseitisnumericallymuchlessstablethancomputing text) is
Γ[G]=TrlogG Trlog((G−1 G−1)G)+EH[ρ]+Exc[ρ]+ΦDMFT[G ] ΦDC[ρ ]+E . (12)
− 0 − loc − loc nuc−nuc
First, we extremize it (δΓ[G]/δG = 0) to obtain the where
Dyson equation
d∆
E =Tr((∆+ε ω )G )
imp imp n imp
G−1 G−1+V +V +Σ V =0. (13) − dωn
− 0 H xc DMFT − dc 1
+ Tr(Σ G ) TS (17)
A note is in order here. We assumed δP/δG = 0, which 2 imp imp − imp
holdswhenevertheprojectordoesnotdependontheself-
Hence
consistent charge density. To ensure this property, we
uthseedrafdoiralthweavloecfaulnizcetdionorubi(tra)lsis|φth(cid:105)e=soluulr(trio)Ynlmof(trˆh),ewschaelarer F = 21Tr(ΣimpGimp)−ΦDC[ρloc]−TSimp
l
relativistic Dirac equation on the LDA charge density +Enuc−nuc Tr((VH +Vxc)ρ)+EH[ρ]+Exc[ρ]
−
(rather than self-consistent charge density). Note also +Trlog(G) Trlog(G )+Tr(V ρ )+µN
loc dc loc
−
that the use of the self-consistently determined Wannier d∆
functions (which depend on self-consistent charge), as is +Tr((∆+εimp ωn )Gimp) (18)
− dω
n
commonly used in most of the LDA+DMFT implemen-
tations[24,29,30],leadstonon-stationaryLDA+DMFT Again using the identity Gimp =Gloc and ρimp =ρloc as
solution, and non-stationary free energies. well as the definition of E2 (Eq. 10) and E3 (Eq. 11) we
We next insert the expression G−1 G−1 into Eq. 12 obtain
− 0
to obtain expression for free energy
F = Trlog(G)+µN +E
2
F =E Tr((V +V )ρ)+EH[ρ]+Exc[ρ] d∆
nuc−nuc− H xc + Tr((∆ ωn +εimp+Vdc)Gloc)
+TrlogG Tr(ΣDMFTG)+ΦDMFT[Gloc] − dωn
−
+Tr(Vdcρloc) ΦDC[ρloc]+µN(14) − Trlog(Gloc)+E3−TSimp (19)
−
This equation is implemented in our DFT+DMFT code.
TheimpurityfreeenergyF containsΦDMFT[G ]in
imp imp Similarly than in the implementation of the total energy
the following way
Eq.8, wecomputeE andTS usingimpurityquanti-
3 imp
F =TrlogG Tr(Σ G )+ΦDMFT[G (].15) ties, while the rest of the terms are computed using the
imp imp imp imp imp
− Green’s function of the solid. [Alternatively, we could
In DMFT, G = G and Σ = Σ , hence we alsocomputethelasttworowsusingimpurityquantities,
loc imp DMFT imp
can write andthefirstrowusingthesolidGreen’sfunction]. Inthis
wayweensurethatF andE aresplitinthesamewaybe-
F =Enuc−nuc Tr((VH +Vxc)ρ)+EH[ρ]+Exc[ρ] tween the ”impurity” and the ”lattice” quantities, hence
−
+Trlog(G) Trlog(G )+F theysharealmostidenticalMonteCarlonoise. However,
loc imp
−
+Tr(V ρ ) ΦDC[ρ ]+µN.(16) when comparing E(V) and F(V) at two different vol-
dc loc loc
− umes,F(V)convergesfasterthatE(V)withthenumber
This equation appears as Eq.4 in the main text. of LDA and/or DMFT iterations. Moreover, F(V) is
Next we split free energy of the impurity into the en- veryrobustwithrespecttosmallchangesinprojectoror
ergy and the entropy term double-counting, while E is more sensitive.
Notice that F + TS can be evaluated at each
imp
F =E TS , LDA+DMFT iteration, just like the total energy above.
imp imp imp
−
9
To add TS at low temperatures, we however need pute is Trlog(G), which requires eigenvalues (but not
imp
a few extra impurity runs. The method of computing eigenvectors) of the LDA+DMFT eigenvalue problem.
TS is explained in the main text, and requires the We first diagonalize
imp
impurity energy at a few temperatures. An alternative
to this approach is to compute TSimp from so called ( 2+V +V +V +Σ(iω ) V )ψ =
”flat-histogram sampling method” [54], which is also −∇ ext H xc n − dc i,k,ωn
=ε ψ .(20)
done as postprocessing on self-consistent LDA+DMFT i,k,ωn i,k,ωn
hybridization ∆.
Perhaps, the most challenging term in Eq. 19 to com- and then evaluate
(cid:88) (cid:88)
Trlog(G)+µN =T (log(εi,k,ωn −iωn−µ)−log(εi,k,∞−iωn−µ))−T log(1+e−β(εi,k,∞−µ))+µN(21)
iωn,i,k,σ i,k,σ
SrVO3
400
LDA+DMFT
350
LDA
300 PBE
eV]250 PBEsol
m
y[200
g
ner150 V
E exp
100
50
0
48 50 52 54 56 58 60 62 64
Volume[ 3]
FIG. 4: Free energy of LDA+DMFT for SrVO compared
3
with total energy of other standard DFT functionals.
Here it becomes apparent that if Σ(iω ) is frequency in-
n
dependent, the first term in the brackets vanishes, while
the second term gives (at T =0) the sum of eigenvalues
FIG. 5: Free energy of LDA+DMFT for FeO compared
(cid:88) with total energy of other standard DFT functionals. Upper
Trlog(G)+µN U=0 θ(ε <µ)ε ,
→ → i,k i,k (lower) panel shows non-magnetic (antiferromagnetic) DFT
i,k,σ calculation. LDA+DMFT results are obtained at 300K in
paramagnetic state.
the well known DFT contribution to the total energy.
estimates it for 1.5%, and PBE overestimates for 0.7%.
COMPARISON WITH STANDARD
FUNCTIONALS Hence predictions of standard functionals in the case of
SrVO are quite in line with standard performance in
3
weakly correlated solids. Perhaps, this is not very sur-
Here we compare total energy of LDA, PBE [55],
prising given that SrVO is a metallic moderately corre-
and PBEsol [56] functionals with the free energy of 3
lated system.
LDA+DMFT.
In most weakly correlated solids, LDA underestimates In FeO (Fig. 5), all standard functionals severally un-
lattice constants on average for 1.6%, while PBE [55] derestimate volume in the paramagnetic state. For ex-
overestimatesthemforapproximately1%.[2]PBEsol[56] ample the lattice constants with LDA, PBEsol and PBE
was designed to predict most accurate volumes in solids, are 7.7%, 6.5% and 5.1% too small, far outside the stan-
and it typically falls in-between LDA and PBE. dard performance of these functionals in weakly corre-
In Fig. 4 we compare LDA+DMFT free energy in lated solids.
SrVO withthetotalenergycomputedbyotherfunction- The predictions are improved when the AFM long
3
als. Both LDA+DMFT and PBEsol underestimate lat- range order is allowed. LDA and PBEsol still underes-
tice constant for approximately 0.6%, while LDA under- timate lattice constant for 3.6%, and 2.3% respectively.
10
TheequilibriumvolumeinCeriumisstronglytemper-
Cerium
100 aturedependent,andisapproximately34˚A3atzeropres-
80 sure and 400K, while it changes to approximately 28˚A3
in the α phase at low temperature. The LDA+DMFT
Energy[meV] 4600 LPPLDDBBEEAAs+oDlMFT i(rsaeslsruoelmatsdeyawrahetac1toGmunPpduae)tretedhsteaimtag4ar0te0eedKm(,e1hn.te5n%wci)et,habteuxptp=uenr0idmteherneptvrieosslcusuomnree-
siderably improved.
20
The DFT results should be compared to T = 0 ex-
020 22 24 26 28 30 32 34 36 perimental volume of 28˚A3. All functionals underesti-
Volume[ 3]
mate the lattice constant, LDA for 6%, PBEsol for 5%
and PBE for 1.8%. Clearly electronic correlations are
FIG. 6: Free energy of LDA+DMFT for Cerium compared
with total energy of other DFT functionals. LDA+DMFT very important even in the α phase at low temperature,
results are obtained at 400K. asstandardDFTfunctionalssubstantiallyunderestimate
the volume.
On the other hand PBE is this time quite close to the
experiment (underestimates for 0.7%). In comparison
LDA+DMFT underestimates it for only 0.16%. It is SCREENED COULOMB REPULSION OF
quite clear that the excellent prediction of AFM-PBE YUKAWA FORM
here is merely a coincidence, as normally PBE overesti-
mates the volume. It is noted above that we used the Yukawa represen-
Finally,weplotresultsforCeriuminFig.6. Theresult tation of the screened Coulomb interaction, in which
of LDA+DMFT is very different from those of any other there is unique relationship between the Hubbard U and
functional, asitclearlycontainsthenontrivialsoftmode Hund’s coupling J. If U is specified, J is uniquely de-
for the α-γ transition. No other functional shows any termined. To show this we derive the matrix elements of
hint of such transition. screenedCoulombinteractioninourDMFTorbitalbasis
(cid:90) (cid:90) (cid:18)u (r)(cid:19)2(cid:18)u (r(cid:48))(cid:19)2 e−λ|r−r(cid:48)|
U = d3r d3r(cid:48) l l Y∗ (ˆr)Y (ˆr)Y∗ (ˆr(cid:48))Y (ˆr(cid:48)) (22)
m1m2m3m4 r r(cid:48) lm1 lm4 lm2 lm3 r r(cid:48)
| − |
There exist a well known expansion of Yukawa interaction in terms of spheric harmonics Y , which reads
km
e−λ|r−r(cid:48)| =4π(cid:88)Ik+1/2(r<)Kk+1/2(r>)(cid:88)Y∗ (ˆr)Y (ˆr(cid:48)) (23)
r r(cid:48) √r< r> km km
| − | k m
Here r = min(r,r(cid:48)), r = max(r,r(cid:48)), I and K are modified Bessel function of the first and second kind. Inserting
< >
this expression into Eq. 22, we get
(cid:88) 4π
U = Y Y Y Y Y∗ Y
m1m2m3m4 2k+1(cid:104) lm1| km1−m4| lm4(cid:105)(cid:104) lm2| km3−m2| lm3(cid:105)
k
(cid:90) ∞ (cid:90) ∞ I (λr )K (λr )
(2k+1) dr dr(cid:48)u2(r)u2(r(cid:48)) k+1/2 < k+1/2 > . (24)
× 0 0 l l √r< r>
Hence, the screened Coulomb interaction has the Slater form with the Slater integrals being
(cid:90) ∞ (cid:90) ∞ I (λr )K (λr )
Fk =(2k+1) dr dr(cid:48)u2(r)u2(r(cid:48)) k+1/2 < k+1/2 > . (25)
0 0 l l √r< r>
Thisisaproductoftwoone-dimensionalintegralsandis It is clear from Eq. 25 that λ uniquely determines
very easy to efficiently implement. all Fk’s, and furthermore even one Slater integral (F0)