Table Of ContentExploring approximations to the GW self-energy ionic gradients
C. Faber1,2,3, P. Boulanger1,2, C. Attaccalite1,2, E. Cannuccia4,5, I. Duchemin2,3, T. Deutsch2,3, X. Blase1,2
1CNRS, Institut N´eel, 38042 Grenoble, France,
2Universit´e Grenoble Alpes, 38042 Grenoble, France,
3INAC, SP2M/L sim, CEA cedex 09,
38054 Grenoble, France.
4Institut Laue Langevin,
BP 156, 38042 Grenoble, France,
5Laboratoire de Physique des Interactions Ioniques et Mol´eculaires,
Aix-Marseille Universit´e/CNRS,
Campus de Saint J´eroˆme, 13397 Marseille, France.
5
(Dated: March 12, 2015)
1
0 Theaccuracyofthemany-bodyperturbationtheoryGW formalismtocalculateelectron-phonon
2 coupling matrix elements has been recently demonstrated in the case of a few important systems.
n However, the related computational costs are high and thus represent strong limitations to its
a widespread application. In the present study, we explore two less demanding alternatives for the
J calculation of electron-phonon coupling matrix elements on the many-body perturbation theory
8 level. Namely, we test the accuracy of the static Coulomb-hole plus screened-exchange (COHSEX)
2 approximation and further of the constant screening approach, where variations of the screened
CoulombpotentialW uponsmallchangesoftheatomicpositionsalongthevibrationaleigenmodes
] areneglected. Wefindthislatterapproximationtobethemostreliable,whereasthestaticCOHSEX
h
ansatz leads to substantial errors. Our conclusions are validated in a few paradigmatic cases:
p
diamond, graphene and the C fullerene. These findings open the way for combining the present
- 60
p many-body perturbation approach with efficient linear-response theories.
m
PACSnumbers: 71.15.Qe,71.38.-k
o
c
.
s I. INTRODUCTION are not yet available within the framework of MBPT.
c Existing GW calculations of EPC matrix elements have
i
s Electron-phononcoupling(EPC)occupiesaprominent beenthereforebasedonthefrozen-phononapproach,ne-
y cessitating a stepwise displacement of the atoms along
roleinvariousfieldsofcondensedmatterphysics,includ-
h
thephononmodeswithanexplicitevaluationoftheelec-
p ing phonon-mediated superconductivity, photoemission
tronic structure on the GW level for each step. This de-
[ bandgaprenormalization,currentcarriersinelasticscat-
manding approach cannot be reconciled with very large
teringorthelifetimeofhotelectrons. Concerningthecal-
1 unit cells and is thus only feasible for zone-center or
culationofEPCmatrixelementsontheabinitiolevel,up
v
tonow, mainlydensityfunctionaltheory(DFT)1 andits zone-boundaryphononmodes. Goingbeyondthefrozen-
8
5 perturbative linear response extensions (DFPT)2,3 have phonon approach, in order to access EPC matrix ele-
ments on the GW level for a large number of electron
0 beenapplied,providingimportantinformationatthemi-
7 croscopic level. and phonon wave vectors, remains a considerable chal-
0 lengetotheab initio community. Significantworkisstill
Recently,severalstudiesquestionedtheaccuracyofthe
. neededtocomeupwithaschemethatmayallowtogen-
1 DFT-based EPC matrix elements when obtained with
eralizethepioneeringabove-mentionedstudiestoamuch
0 (semi)local functionals such as LDA or PBE.4 By way
5 larger set of systems.
of example, the electron-phonon coupling involving spe-
1
cific electronic and phonon modes in graphene5,6 and In the present study, we explore the merits
:
v graphite,7 the value of the electron-phonon coupling po- of the static coulomb-hole plus screened exchange
Xi tential to states at the Femi level in the electron-doped (COHSEX)14,15,19–27 approximation to the GW self-
fullerene,8,9superconductingbismuthatesandtransition- energy for the calculation of the electron-phonon cou-
r
a metalchloronitrides,10 ortherenormalizationofthepho- pling matrix elements in the isolated fullerene molecule,
toemission band structure of pentacene11 and diamond in diamond and in graphene. Our definitions of the rele-
crystals,12 were shown to be affected by a significant un- vantEPCmatrixelementsfollowpreviousstudiesrelated
derestimation of the EPC matrix elements when calcu- to the superconducting transition in the fullerides,8,9
lated within DFT and (semi)local functionals. the zero-point motion renormalization of the gap in
As a cure to such problems, many-body perturba- diamond12 andthegapopeningthroughelectron-phonon
tiontheory(MBPT)techniqueswithintheso-calledGW coupling in graphene,5 respectively. For the same set of
approximation13–18 showed a clear improvement when systems, we also explore the accuracy of the constant-
compared to available experimental data.5,7,9–12 Unfor- screening approximation, namely the assumption that
tunately, the kind of theories that are available within thescreenedCoulombpotentialcanbeconsidered,tolin-
DFT,andinparticularthepowerfulDFPTformalism,2,3 earorder,asaconstantuponsmallchangesoftheatomic
2
positions. WhereasCOHSEXleadstonon-negligibledis-
crepancies,thislatterapproximationwillbeshowntobe
occ
robust and accurate. The present results offer impor- ΣSEX(r,r(cid:48);E) = (cid:88)W(r,r(cid:48);E ε )φ (r)φ∗(r(cid:48))
tantperspectivesinmakingtheGW formalismamenable − − n n n
n
to the study of condensed matter phenomena involving (cid:88) (cid:90) ∞ B(r,r(cid:48);ω(cid:48))
electron-phonon coupling. ΣCOH(r,r(cid:48);E) = φ (r)φ∗(r(cid:48)) dω(cid:48) ,
n n P E ε ω(cid:48)
n 0 − n−
II. METHODOLOGY where indicates the principal value. ΣSEX, which in-
P
volves a summation over the occupied states only, orig-
inates from the poles of G. It is called the screened ex-
Webrieflyintroducethemany-bodyperturbationthe-
change interaction in analogy to the bare exchange term
ory Green’s function formalism, providing a solid frame-
that can be obtained by replacing W with the energy-
work for the calculation of quasiparticle energies E. In
independentbareCoulombpotentialv. ΣCOH originates
such an approach, the one-body quasiparticle eigenvalue
from the poles of W and represents the Coulomb-hole
equation reads:
contribution, since it can be shown to be related to the
interaction of an electron with its related adiabatically
(cid:18) 2 (cid:19) built correlation hole.
−∇ +Vion(r) φ(r)+VH(r)φ(r)
2 The static approximation to the exact COHSEX de-
(cid:90) compositionassumesthat(E εn) 0forall(n),leading
+ dr(cid:48) Σ(r,r(cid:48);E)φ(r(cid:48))=Eφ(r), to simplified static screened e−xchan(cid:39)ge and Coulomb-hole
expressions, namely:
where Vion and VH are the ionic and Hartree poten-
tial, respectively. The self-energy Σ(r,r(cid:48);E) replaces occ
(cid:88)
the well-known exchange-correlation potential of den- ΣSEX(r,r(cid:48);0) = W(r,r(cid:48);0)φ (r)φ∗(r(cid:48))
static − n n
sity functional theory or the exchange operator in the n
Hartree-Fock formalism. In general, it is non-local, 1(cid:88)
energy-dependent and non-Hermitian. Derived within ΣCstOatHic(r,r(cid:48);0) = 2 φn(r)φ∗n(r(cid:48))W(cid:102)(r,r(cid:48);0)
Schwinger’s functional derivative approach to perturba- n
1
tion theory,13 the GW approximation to the self-energy = δ(r r(cid:48))W(cid:102)(r,r(cid:48);0).
leads to: 2 −
Here, W(cid:102) =(W v)isthedifferencebetweenthescreened
i (cid:90) and bare Coulo−mb potential. As a result, besides being
ΣGW(r,r(cid:48);E)= dωeiδωW(r,r(cid:48);ω)G(r,r(cid:48);E+ω),
2π a static approximation, this Coulomb-hole term is also
local in space. Such a static COHSEX approximation
whereGandW arethetime-orderedone-particleGreen’s (labeled COHSEX here below) was shown to yield too
function and the dynamically screened Coulomb poten- large gaps in the case of semiconductors.15 By way of
tial, respectively, and δ = 0+ a small positive infinitesi- example, in the present case of the C molecule, the
60
mal. COHSEX gap is found to be 5.3 eV, i.e. 0.4 eV larger
TheCoulomb-holeplusscreened-exchange(COHSEX) than the 4.9 eV experimental gap.67 Nevertheless, this
representationoftheGW self-energywasoriginallyintro- has to be∼compared to the starting 1.6 eV DFT-LDA
ducedusingatimerepresentationofG, W andΣ.14,20,21 Kohn-Sham gap which is dramatically too small.
We follow Hybertsen and Louie15 by using the following While it cannot be claimed that the static COHSEX
spectral representations of G and W: approach is a good approximation to absolute quasipar-
ticle energies, we emphasize that we are interested in
quasiparticle energy differences upon small (infinitesi-
G(r,r(cid:48);E+ω) = (cid:88) φn(r)φ∗n(r(cid:48)) mal) atomic lattice motions. The main assumption on
E+ω ε iδsgn(ε µ)
n − n− n− which we rely to calculate the electron-phonon coupling
(cid:90) ∞ 2ω(cid:48)B(r,r(cid:48);ω(cid:48)) within the COHSEX approximation is that the varia-
W(r,r(cid:48);ω) = v(r,r(cid:48))+ dω(cid:48) ,
ω2 (ω(cid:48) iδ)2 tions of the dynamical contribution to the self-energy
0 − − can be neglected. This can be a priori rationalized by
where G has been written in terms of one-body eigen- emphasizing that dynamical interactions are driven by
states φ and eigenenergies ε , typically starting DFT the plasmons dynamics, collective excitations less sensi-
n n
Kohn-Sham or Hartree-Fock solutions, and W in terms tive to small atomic displacements than single-particle
of its spectral function B. These expressions allow to excitation energies and wave functions. It remains, how-
obtain the pole structure of both G and W. From the ever, that besides the static approximation, the spatially
residue theorem, one then rapidly obtains an exact de- local character of the static COH term is at odds with
composition of Σ: the nonlocality of the full GW self-energy.
3
Inordertofurtherjustifythefollowingresultsconcern- the superconducting transition temperature through the
ing the constant screening approach, namely the second dimensionless parameter λ = N(E )Vep, where N(E )
F F
approximation we explore in this study, it is instructive is the density of states at the Fermi level. Namely, Vep
to consider the GW plus Bethe-Salpeter formalism.28–30 reads in the so-called molecular limit:44,45
Thelatterisamany-bodyperturbationtheoryapproach
concerned with describing the linear response of a sys-
tteiomnwanitdhtrheuspsetchtetoBaetthiem-Se-adlpepeteenrdeeqnutaetxitoenrn(aBlSpEe)rtiusrtbhae- Vep = 9M1 (cid:88)ω12 (cid:88)3 (cid:12)(cid:12)(cid:12)(cid:12)∂∂(cid:15)umν(cid:12)(cid:12)(cid:12)(cid:12)2. (1)
MBPT analogue to time-dependent DFT. At the heart ν ν m=1
of the GW/BSE approach lies the variation (∂GW/∂λ),
Here, M is the mass of the carbon ions, ω is the fre-
wherethe“perturbation”λistheone-bodyGreen’sfunc- ν
quency of the vibrational mode with index (ν) and uν is
tion G. The most common approximation that has been
the vibrational polarization vector. As already discussed
shown to be remarkably accurate31–33 is to replace the
before, the electronic states with energy ((cid:15) ) are lim-
GW self-energy by its static COHSEX approximation, m
ited to the 3-fold degenerate LUMO level. From group
and to consider further that (∂W/∂λ) = 0. It is such a
theory analysis it follows that only the H and A vi-
simplifiedschemeweaimtoexploreinthepresentstudy, g g
brational modes can couple to these states. A schematic
differing in the fact that the perturbation λ is now in-
representation of the t level splitting as a function of
duced by a vibrational distortion of the system. 1u
the deformation amplitude along a H mode is provided
g
in Fig. 1(a). The EPC matrix elements are consequently
relatedtotheslopes(∂(cid:15) /∂uν),showingastrongdepen-
m
III. TECHNICAL DETAILS denceontheformalismadoptedtocalculatetheenergies
(cid:15) . These energy derivatives are calculated within the
m
The many-body GW and COHSEX calculations for frozen phonon approach, using a symmetric five points
theisolatedC fullereneareperformedusingtheFiesta finite-difference formula.
60
code, an implementation of the GW formalism within Such an effective electron-phonon potential was calcu-
a Gaussian basis.34–36 We start from DFT-LDA eigen- latedforC inRefs.8,9,bothwithinDFTusingdiverse
60
statescalculatedwiththeSiestapackage37 andatriple- (semi)local and hybrid exchange-correlation functionals
zeta plus polarization (TZP) basis38 for the description and GW. While DFT-LDA values were shown to signif-
of the valence orbitals combined with standard norm- icantly underestimate the coupling energy, DFT calcula-
conserving pseudopotentials.39 We use the resolution of tionsperformedwithglobalhybridscontaining20%-30%
the identity technique (RI-SVS) with an even-tempered ofexactexchangeandtheGW approachwerefound8,9 to
auxiliary basis set formed of four Gaussians per angular compare favorably to gas phase experimental data.46,47
channel (up to d character) with exponents between 0.2 As a drawback of DFT calculations with hybrids, it was
and 3.2 Bohr−2. With such a basis, our G W @LDA evidenced thatthe resulting Vep potential wouldquickly
0 0
ionization potential and HOMO-LUMO gap are found increasewiththeamountofexactexchange, yieldingthe
to be 7.33 eV and 4.39 eV (B3LYP geometry), respec- usual question of the proper choice of the functional pa-
tively, starting from a 1.66 eV DFT-LDA gap. This is in rameters for a given system.
closeagreementwiththe7.45eVand4.40eVplanewave In the present study, we are concerned with explor-
G W @LDA values of Ref. 40 (LDA geometry). The ing simplified GW schemes, i.e approximations reducing
0 0
G W @LDA gap remains smaller than the 4.9 eV ex- the needed computational effort. As such, quasiparti-
0 0
∼
perimental value, an issue that will be addressed by the cle energies are evaluated at a single-shot COHSEX and
partiallyself-consistentschemeusedinthepresentstudy, a partially self-consistent evCOHSEX level, where the
as discussed below. quasiparticle energies are reinjected self-consistently in
Following Ref. 8, we use the relaxed structure and the construction of G and W, while keeping the start-
phonon eigenmodes generated within a DFT-B3LYP ap- ing Kohn-Sham LDA wave functions unchanged. Our
proach using a 6-311G(d) basis. This was shown to pro- reference point are partially self-consistent evGW calcu-
vide vibrational frequencies in excellent agreement with lations with self-consistency on the eigenvalues as pre-
Raman measurements.8,41 Concerning the calculation of sented in Ref. 9 for the calculation of the EPC matrix
the EPC matrix elements, our present study follows elements in C . Such a simplified approach to full self-
60
the approach described in Refs. 8,9, where the electron- consistency,justifiedbythedramaticallytoosmallstart-
phonon coupling strength associated with the (electron- ing point DFT-LDA gap for small isolated molecules,
doped) lowest unoccupied 3-fold degenerate t molecu- has been shown in the case of gas phase organic systems
1u
lar orbital (LUMO) has been explored within DFT and to produce more accurate ionization potentials and elec-
GW. The specific choice of the LUMO level is dictated tronicaffinities,34,35,48–51 togetherwithimprovedoptical
by the physics of the phonon-mediated superconducting excitation energies at the GW/BSE level,36,52–54 when
transition in the fullerides.42,43 We adopt the standard LDAorPBEKohn-Shameigenstatesarechosenasstart-
definition of the effective electron-phonon coupling po- ing points. Consistently, when starting from DFT-LDA,
tential Vep that enters e.g. the McMillan formula for EPC matrix elements in C were also found to be in
60
4
FIG.1: Schematicrepresentationofchangesintheelectronicstructureunderphonondistortion: a)inC ,the3-foldLUMO
60
levels (blue) split under a distortion (red) along a H phonon mode. b) In graphene, an energy gap ∆E at the Dirac point is
g
opened(red)whendistortingtheequilibriumstructure(blue)alongtheΓ−A(cid:48) orK−E phononmodes. c)Indiamond,the
1 2g
lowestlyingconductionbandsatΓ/X (blue)arerenormalisedbytheircouplingtoazone-boundaryphononmode. Forclarity,
only the states showing the largest shift are plotted (red dot, see text).
closer agreement with experiment within evGW as com- tive (non self-consistent) GW approach to obtain quasi-
pared to G W .9 particle energies.
0 0
For the periodic diamond and graphene systems, we For diamond, following Antonius and coworkers,12 we
use the Yambo code which implements many-body per- haveusedarelaxedDFT-LDAgeometry(latticeparame-
turbation theory within a plane wave formalism.55 The tera=3.591˚A)withaplanewavecutoffenergyof80Ry
startingDFTKohn-Shameigenstatesaregeneratedwith and norm-conserving pseudopotentials. All calculations
the Abinit code.56 In the case of graphene, the running were performed with a 8x4x4 uniform Monkhorst-Pack
parameters are identical to those used in a previous GW grid to describe the supercell necessary for the frozen-
study of the electron-phonon coupling in this material,5 phonon simulation of the (q = X) phonon. We have
namely we use a = 2.46 ˚A as lattice parameter, 15 a.u. used 160 bands and a cutoff of 8 Ha for the construc-
distance between the layers,57 36 36 1 k-points for tion of the dielectric function and the GW self-energy.
× ×
the GW calculations in the Brillouin zone (BZ) corre- Following Refs. 12 and 58 where the diamond band-gap
sponding to the primitive (2-atoms) cell and the nearest renormalization by the zero-point motion was studied,
equivalent k-grid in the BZ corresponding to the super- we define the relevant EPC matrix elements as the sec-
cell needed to describe zone-boundary (q=K) phonons ond derivative of the electronic eigenvalues with respect
inthefrozen-phononapproach. Theenergycutoffonthe to the phonon displacements at equilibrium:
plane wave basis was set to 60 Ry. All conduction bands
within an energy range of 45 eV above the Fermi level
wereusedtobuildGandW,andacutoffof4Haforthe ∂(cid:15) ¯h ∂2
mk = (cid:15) , (3)
constructionofthedielectricfunction. TheGodby-Needs ∂n 2Mω ∂u2 mk
νq νq νq
plasmon pole model was employed for the dynamical di-
electric constant entering in W.
where n is the occupation of the phonon state (ν) of
νq
As schematically described in Fig. 1(b), we focus on
wave vector q. The EPC matrix elements are now re-
the splitting ∆E of the degenerate occupied and unoc-
lated to the temperature dependent renormalization of
cupied levels at the Dirac point caused by the coupling
the electronic bands. We evaluate the second-derivative
to the qν =Γ-E optical mode at zone center together
2g of the eigenvalues using a symmetric five-point finite-
with the zone-boundary qν = K-A(cid:48) phonon, for which difference formula with displacements of 0.01, 0.02 and
1
the strongest renormalization of the coupling constant 0.03 ˚A. In Ref. 59, it was already shown that stan-
have been observed upon replacing the DFT-LDA ap- dardDFTandDFPTmethodswith(semi)localfunction-
proach by the G0W0 formalism.5,6 The EPC matrix ele- als strongly underestimate (by ca. 30%) the zero-point
mentsofinterestarethushererelatedtothederivatives: renormalization of the direct gap of diamond compared
(∂∆E/∂uqν), namely (see Ref. 5): to experiment, while Ref. 12 demonstrated that G0W0
calculationsrestoredthegoodagreementbetweentheory
andexperiment. Wefocusherebelowontherenormaliza-
1 (cid:18)∂∆E(cid:19)2
D2 = , (2) tion of the lowest conduction band edge at (mk = X1c)
(cid:104) qν(cid:105) Aqν ∂uνq and (mk = Γ15c) by the coupling to the (νq = X4)
phonon modes at the zone-boundary (see Fig. 1c). It
expressed in (eV/˚A)2, where Aqν equals 16 and 8 for is for such a phonon wave vector that the EPC matrix
qν = Γ-E and qν = K-A(cid:48), respectively. We per- elements were found to be large and that also the GW
2g 1
formed“single-shot”G W andCOHSEXcalculationson correction was shown to be the most significant.12 Fol-
0 0
top of DFT-LDA results, namely the standard perturba- lowing Ref. 12, the provided EPC matrix elements are
5
LDA G W G W (W) COHSEX COHSEX(W)
0 0 0 0
Graphene
Γ-E 44.3 68.1 (65.2) 72.0 81.5 83.2
2g
(cid:48)
K-A 88.6 201(187) 207 260 256
1
Diamond
Γ 1.281 0.611 0.657 1.311 1.445
15c
X -2.009 -1.031 -1.066 -2.124 -2.273
1c
TABLE I: Calculated EPC-related quantities (see text) for graphene and diamond. In the case of graphene, these correspond
to the splitting of the degenerate highest valence and lowest conduction band at the Dirac point under the influence of a
(qν = Γ-E ) and (qν = K-A(cid:48)) phonon mode and are expressed in (eV/˚A)2. Number in parenthesis are obtained with a
2g 1
plasmonpoleenergyof7eV,closetothegrapheneπ-plasmonenergy,insteadofthe27eVoriginallyusedinRef.5(seetext).
Fordiamond,thesearetheresultsofequation(3)fortheX andΓ conductionstatesandfortheX zone-boundaryphonon
1c 15c 4
mode only (values in eV). This corresponds to twice the zero-point renormalization arising from this given phonon mode (see
Ref. 12).
averaged over the degenerate electronic manifold consid- tronic state cannot be made. Our findings are, however,
ered (namely the 3 degenerate levels at zone-center and very consistent with what is observed for the Γ state,
15c
the 2 degenerate levels at the X zone boundary). namelyadramaticdecreaseofthecouplingstrengthwith
1c
the X phonon mode which, again, is the dominant cou-
4
pling mode. The coupling energy is indeed found to be
IV. RESULTS reduced by 52% and 49% from LDA to GW for the Γ
15c
andX states,respectively. SincetheX stateiscloser
1c 1c
A. Reference GW calculations to the true conduction band minimum, this is a strong
indication that zero-point motion renormalization of the
indirect gap will certainly also be strongly affected by
We first reproduce for validation and for reference the
the GW correction. The full study of such an effect is
previously published GW calculations related to EPC
beyond the scope of the present paper.
matrix elements in graphene, diamond and C . Our re-
60
For C , the total evGW coupling potential is within
sults are compiled in Table I for diamond and graphene 60
8% of that found by Ref. 9 as a result of the larger TZP
andinTableIIforC ,respectively. OurLDAandG W
60 0 0 basis used here.38 Our evGW Vep potential is found to
values for graphene are close to those found in Ref. 5,
namely e.g. 201 eV/˚A (G W value, present study) for increase by 44% with respect to the corresponding LDA
0 0
the largest matrix element with the (K-A(cid:48)) phonon, to value, to be compared with the 48% increase obtained
1 in Ref. 9 with a smaller basis. This indicates the good
be compared to 193 eV/˚A in the previous study.57 The
convergenceoftheGW correctiontotheLDAvalue. We
difference can be explained by an increased five points,
note that while DFT and GW calculations are known
instead of only two, finite-difference formula, and by the
toconvergeverydifferentlywithrespecttobasissize, we
increaseofthedielectricmatrixsizeenergycutoff,from2
compareinthepresentstudyapproacheswhicharemuch
to4Hacutoff. WealsoverifiedthatchangingtheGodby-
moresimilar,namely“standardandapproximated”GW
Needs plasmon model input finite frequency, from the
calculations, indicating certainly an even better conver-
(default) 27 eV value in Ref. 5 to a 7 eV value closer to
gence of the differences we are interested in.
the π-plasmon resonance in graphene, does not change
significantly the calculated EPC strength (see numbers
in parenthesis in the Table I). Such differences are neg-
ligible with respect to the more than 100% increase as B. The COHSEX approximation
compared to the LDA value.
For diamond, the EPC contribution associated with We now explore approximations to the full GW cal-
the mk = Γ conduction states are in good agree- culations performed in Refs. 5,9,12 and reproduced here
15c
ment with the LDA and G W results of Antonius and above. Westartbythecaseofdiamond. Ascomparedto
0 0
coworkers,12 i.e. we find a 11 meV difference with re- the G W calculations, the static COHSEX approxima-
0 0
spect to their EPC in the case of LDA, while we find a tion is shown to induce errors larger than 100% for the
difference of 8 meV for the G W results.68 These small relevantelectron-phononcouplingenergytotheΓ and
0 0 15c
variations can be ascribed to the different pseudopoten- X states. In both cases, the static COHSEX approach
1c
tials (Ref. 59 showed that pseudopotentials can lead to strongly overestimates the electron-phonon coupling. A
errors up to 50 meV), and to the different convergence representation of the errors as compared to the G W
0 0
parameters. calculations is provided in Fig. 2.
SinceRef.12focusedontherenormalizationofthedi- The case of graphene is of specific interest since, in
rect band gap at the zone-center, comparison with the greatcontrasttodiamondandC ,thefundamentalgap
60
present results for the phonon coupling to the X elec- reduces to zero in the equilibrium geometry and opens
1c
6
Theory Experiment
Mode LDA evGW COHSEX evCOHSEX evCOHSEX(W)
Fullerene
H (1) 4.7 5.85 4.99 6.3 6.3
g
H (2) 9.3 10.1 9.99 10.3 11.9
g
H (3) 8.8 12.95 11.96 13.7 15.3
g
H (4) 4.2 5.5 4.9 5.4 5.7
g
H (5) 3.98 4.9 4.5 5.3 5.5
g
H (6) 1.8 1.99 2.05 2.2 2.4
g
H (7) 15.8 25.4 23.3 27.9 26.5
g
H (8) 13.1 18.3 17.1 19.8 19.5
g
A (1) 1.2 1.9 1.85 1.8 1.2
g
A (2) 7.2 13.98 13.2 15.4 12.2
g
Total A 8.4 15.9 15.0 17.2 13.4
g
Total H 61.7 84.96 78.8 90.9 93.1 96.2b,96.5c
g
Total 70.2 100.9 (108.6a) 93.8 108.1 106.6 106.8b
TABLEII:Mode-resolvedVep couplingpotentialassociatedwiththe3-folddegenerateC LUMOandthecorrespondingVep
ν 60
totals (in meV). For group symmetry reasons, only the H and A vibrational modes couple to this electronic state.
g g
a Ref. 9,
b Ref. 46, Table V,
c Ref. 47.
Diamond Graphene
350
300 Γ15cstate W) X1cstate W) 40 Γ(E2g)mode K(A1’)mode
%err112205050000 LDA W) COHSEX COHSEX( LDA W) COHSEX COHSEX( %err−22000 LDA GW00 W(W)0 HSEX EX(W) LDA GW00 (W) HSEX EX(W)
50 GW00 GW(00 GW00 GW(00 −40 G0 CO COHS GW00 CO COHS
0
60
−
FIG. 2: Relative error of the different approximations with FIG. 3: Relative error of the different approximations with
respect to the G0W0 results for the “thermal renormaliza- respecttotheG0W0 resultsfortheelectron-phononcoupling
tion” EPC energies (Eqn. 3) for the lowest conduction band in graphene.
of diamond.
is observed again to systematically increase the coupling
upon phonon distortions. For the optical phonon mode
strength as compared to evGW.62 This result is consis-
at the zone-center, the COHSEX approximation leads
tent with what was observed for diamond and graphene,
again to an increase as compared to the G W value,
0 0 even though clearly the error appears to vary from one
eventhoughsmallerthaninthediamondcasewitha20%
system to another.
error. Such an error becomes larger for the K-A(cid:48) mode
1
witha29%increaseoftheEPC-relatedgap-openingrate.
In the present graphene case, the error induced by the
COHSEX approach remains smaller than the 34% and C. The constant screening approximation
56% reduction observed upon using DFT-LDA instead
of G W . However, the close to 30% error observed for Following the above-mentioned analogy with the
0 0
the COHSEX calculation of D2 is clearly significant Bethe-Salpeter formalism, we finally test another ap-
(cid:104) K(cid:105)
(see Fig. 3), questioning again the applicability of the proximation, that is the constant screening approach.
static COHSEX formalism for EPC matrix elements cal- Namely, we explicitly calculate the screened Coulomb
culations. potential for the undistorted structure and upon chang-
We now consider the C case, where we compare ing the positions of the ions along the (R ) vibrational
60 qν
evGW and evCOHSEX calculations. A global increase eigenmodes, we keep W frozen to the undistorted value.
of about 7% from evGW to its static approximation is Formally,thisamountstoassumingthat: ∂GW/∂uqν
observed. Besides the H (4) and A (1) modes showing (∂G/∂uqν)W. SuchanapproachwillbelabeledGW(W≈)
g g
verysmallcouplings,thestaticCOHSEXapproximation or COHSEX(W) in the following.
7
Fullerene Comparingthetotalcoupling,evCOHSEX(W)andev-
30
COHSEX agree within 1.5 %.
20 Hgmode
We now explore the constant-screening approximation
10 for periodic diamond. Results are again compiled in Ta-
0 ble I. For diamond, we observe a 3.4% and 7.5% dif-
%err−−2100 evGW COHSEX HSEX(W) cEXP dEXP fcΓeo1rn5ecsntcasetnatbtseecstrw,eeereennsinpGgec0atWipvpe0rl(yoW.xi)mAaatntidotnhGey0iWCelOd0sHfeoSrrErotXrhseloeXfvet1hlc,eatonhrde-
30 v O der of 7.0% and 10.2%, respectively. In summary, while
− A e C
−40 LD ev the COHSEX approximation was failing for diamond,
the constant-W approximation leads again to very rea-
50
− sonnable results, with errors well within 10% at the GW
level. As illustrated in Fig. 2, this error is much smaller
FIG. 4: Relative error of the different approximations with
respect to the evGW results for the H -related Vep of equa- than the error induced by the LDA approximation.
g
tion (1). For convenience, a comparison to the experimental We finally address the case of graphene. Concern-
data is provided. ing the coupling with the zone-center optical mode, the
constant-W approach leads to a 5.7% and 2.1% error
respectively when comparing G W (W) to G W and
0 0 0 0
Inthecaseofthefullerene,wherecalculationsareper- COHSEX(W) to COHSEX. In the case of coupling with
formed within a Gaussian basis, special care must be the zone-corner K-A(cid:48) phonon mode, the discrepancy is
1
takenintheimplementationofsuchaconstant-screening of 3% and 1.5% applying the constant-screening approx-
approximation. Intheauxiliarybasis,orresolutionofthe imation to GW and COHSEX approaches, respectively.
identityapproach,thebareandscreenedCoulombpoten- Notice that the results for the constant-screening of the
tialsareexpressedintermsofanatom-centeredauxiliary K-A(cid:48) phonon are affected by a numerical error of about
1
basis with the following relations: 4%, therefore the discrepancy could be slightly larger.
Anywaythisdifferenceismuchsmallerthantheerrorin-
(cid:90) (cid:90) duced by the standard DFT-LDA calculations, but also
[W]β,β(cid:48) = drdr(cid:48)β(r)W(r,r(cid:48))β(cid:48)(r(cid:48)), (4) much smaller than the error induced by the static COH-
SEX formalism as compared to the reference GW calcu-
W(r,r(cid:48))=(cid:88)β(r)(cid:0)S−1[W]S−1(cid:1) β(cid:48)(r(cid:48)), (5) lation.
β,β(cid:48)
β,β(cid:48)
where S is the overlap matrix in the auxiliary basis. Us-
V. DISCUSSION
ingnowthenotation: W andβ forthescreenedCoulomb
potential and the auxiliary basis for the slightly dis-
torted system, the assumption: W(r,r(cid:48)) W(r,r(cid:48)) Clearly, among the two approximations explored here
leads straightforwardly to the condition: (cid:39) above, the constant-screening stands as a much better
approach than the static COHSEX for the calculation
of electron-phonon coupling matrix elements within the
[W] S [W] S , (6) GW formalism. Overall, the largest error induced by
β,β(cid:48) ββ β,β(cid:48) β(cid:48)β(cid:48)
(cid:39)
the static COHSEX approximation is larger than 100%
where S =< β β > is an overlap matrix between the in the diamond case, while it is reduced to 7% in the
ββ
|
auxiliary bases for the perturbed and unperturbed sys- caseoftheconstant-screeningapproximationtofullGW
tems, respectively. calculations.
For C again, we test the constant screening approx- The static COHSEX approximation is known to al-
60
imation at the COHSEX level only, namely comparing ways overestimate band gaps. As proposed in Ref. 8, a
evCOHSEXand evCOHSEX(W)calculationsof theVep too large band gap should lead to an underscreening of
energies. While the plasmon-pole approach used for dia- the electron-phonon interaction and consequently to en-
mond and graphene is based on a fixed finite frequency hancedEPCmatrixelements. Thisisconsistentwiththe
pole value independent of the geometry, the real axis fact that, within Hartree-Fock or with increasing exact
poles contribution to the correlation energy in the con- exchangeinhybridfunctionals,theEPCmatrixelements
tourdeformationapproach(seeRefs.19,34)implemented are found to steeply increase. The Hartree-Fock D2
intheFiestapackagechangesfromonestructuretoan- and D2 coupling constants in graphene were foun(cid:104)dΓto(cid:105)
(cid:104) K(cid:105)
other,leadingtodifficultieswhentryingtoimplementthe be about 5 and 30 times larger, respectively, as com-
constant-W approach within the frozen-phonon scheme. pared to their GW analog,5 while in C increasing the
60
The results of the constant-screening approximation percentage of exact exchange from 20% to 30% in hy-
are compiled in Tables I and II. For C , an excel- brid functionals was found to enhance the Vep energy
60
lentagreementisobtainedwithinevCOHSEX(W)com- by about 15%.8 Such an interpretation matches the ob-
pared to the corresponding evCOHSEX calculations. servation that in graphene and C , the COHSEX EPC
60
8
coupling constants are larger than their corresponding VI. CONCLUSIONS
GW reference. In the case of diamond, where we are
concernedwiththesecond-orderderivatives,theeffectof
assuming a static approximation may be more difficult We have explored two approximations for calculating
to interpret. self-energy gradients of interest for electron-phonon cou-
AnimportantobservationfurtheristhattheCoulomb- pling, namely the two simplifications commonly used in
hole term is a local potential in the static COHSEX ap- the GW/Bethe-Salpeter calculations, that is the static
proximation, washing out spatial local-field effects, a ob- COHSEX and the constant screening approximations.
servation that may potentially allow to understand the We explored these approaches in the case of diamond,
failure of the COHSEX approximation to reproduce ac- graphene and C . Our findings suggest that the COH-
60
curatelytheevolutionofthebandstructurewithphonon SEXapproximationcannotbetrustedtoimproveonthe
deformation. In any case, the large variations of the in- DFT-LDA values as clearly illustrated in particular in
duced errors from diamond to graphene and C60 is still the case of diamond. On the contrary, the constant
to be understood, and we will just stand here on the ob- screening hypothesis, namely assuming that W remains
servationthatthestaticCOHSEXapproximationcannot to first order constant with respect to small ionic dis-
be trusted to improve on DFT-LDA calculations. placements, seems to be reliable, with a discrepancy no
Concerning the neglect of the gradient of W with larger than 10%, even in the difficult case of graphene
respect to the ionic positions, such an approximation where the phonon perturbation dramatically affects the
was already tested by Ismail-Beigi and Louie in the Diracconeandthe(semi)metallicnatureofthegraphene
context of excited state ionic forces within the Bethe- sheet. Even though a deeper understanding of the ori-
Salpeter formalism,63 showing good agreement with ex- gin of the specific difficulties uncountered by the static
plicit finite-difference ”exact” BSE calculations for small COHSEX approximation would allow to better rational-
CO and NH3 molecules. It is commonly assumed that ize the validity and limits of the tested approximations,
the GW/BSE approach is much more resistant to ap- the present results offer promising perspectives to carry
proximations on W as compared to the GW approach on such many-body evaluations of the electron-phonon
for (charged) excitations. This is due to cancellations coupling gradients with much reduced computer cost on
of errors between the electron-electron and electron-hole realisticsystems,includingthestudyofperiodicsystems
interactions. Namely, any error introduced in W is ex- witharbitrarywavevectorperturbation. Furtherstudies
pected to affect excitonic interactions and quasiparticle arehoweverrequiredonalargersetofsystemsandphysi-
gapsinoppositeways. Clearly, thepresentGW studyof calobservablesinordertobetterassesstheinterest,with
thevariationsofagivenquasiparticleenergywithrespect respect to common DFT calculations, of using the GW
toionicpositionscannotbenefitfromsuchcancellationof approach,anditsvariousapproximations,forcalculating
errors. Still, the constant-screening approximation turns self-energy gradients with respect to ionic positions.
to be a reliable approach to save on computational cost.
An important consequence of the present findings is Acknowledgments. P.B. warmly acknowledges G.
that once the screened Coulomb potential W(r,r(cid:48);ω) is Antonius for providing the diamond phonon eigenmodes
built for the equilibrium geometry, the calculation of used in this study and for guiding our comparison with
the variations of the quasiparticle energies with respect Ref.12. C.A.acknowledgesDanieleVarsanoforthecor-
to the perturbation (λ) only requires the evaluation of rectionofsomebugsintheYambocode. C.F.isindebted
the variations of the Green’s function G with respect to to the French CNRS and CEA for Ph.D funding. P.B.
the perturbation. This can be performed within stan- acknowledges a postdoctoral fellowship from the French
dard DFPT techniques, at least in the case of non-self- National Research Agency under Contract No. ANR-
consistentG W calculationswheretheGreen’sfunction 2012-BS04PANELS.Computingtimehasbeenprovided
0 0
assumes an explicit form in function of the input DFT by the “Curie” national GENCI-IDRIS supercomputing
eigenstates. This may invite, for systems such as C , to center under contract No. i2012096655 and a PRACE
60
useG W calculationsstartingfromDFTeigenstatesob- European project under contract No. 2012071258. X.B.
0 0
tainedwithhybridfunctionals,whichhavebeenshownto acknowledges illuminating discussions with Steven G.
be a better starting point for finite-size systems.54,64–66 Louie and G. Antonius.
1 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864-B871 5 M. Lazzeri, C. Attaccalite, L. Wirtz, and F. Mauri, Phys.
(1964). Rev. B 78, 081406 (2008).
2 S.Baroni,S.deGironcoli,A.DalCorso,andP.Giannozzi, 6 D. Basko and I. L. Aleiner, Phys. Rev. B, 77, 041409,
Rev. Mod. Phys. 73, 515 (2001). (2008).
3 X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997). 7 A. Gru¨neis, J. Serrano, A. Bosak, M. Lazzeri, S. L.
4 J.P.Perdew,K.BurkeandM.Ernzerhof,Phys.Rev.Lett. Molodtsov, L. Wirtz, C. Attaccalite, M. Krisch, A. Ru-
77, 3865-3868 (1996). bio, F. Mauri, and T. Pichler, Phys. Rev. B 80, 085423
9
(2009). descriptionoftheeigenstates.Thisbasisislargerthanthe
8 J. Laflamme Janssen, M. Coˆt´e, S.G. Louie, and M.L. Co- double-zeta plus polarization (DZP) basis used in Ref. 9,
hen, Phys. Rev. B 81, 073106 (2010). withtheinclusioninparticularofdiffuses andp orbitals.
9 C. Faber, J. Laflamme Janssen, M. Coˆt´e, E. Runge, and 39 N. Troullier and J.L. Martins, Phys. Rev. B 43, 1993
X. Blase Phys. Rev. B 84, 155104 (2011). (1991).
10 Z.P. Yin, A. Kutepov, and G. Kotliar, Phys. Rev. X 3, 40 T. Anh Pham, H.-V. Nguyen, D. Rocca, and G. Galli,
021011 (2013). Phys. Rev. B 87, 155148 (2013).
11 S.Ciuchi,R.C.Hatch,H.Ho¨chst,C.Faber,X.Blase,and 41 Onehastonotehoweverthattheimpactofthefunctional
S. Fratini, Phys. Rev. Lett. 108, 256401 (2012). on the vibrational eigenmodes used to calculate the EPC
12 G. Antonius, S. Ponc´e, P. Boulanger, M. Coˆt´e, and X. matrixelementswasfoundtobemarginal,themaineffect
Gonze, Phys. Rev. Lett. 112, 215501 (2014). comingfromtheframeworkwithinwhichthequasiparticle
13 P.C. Martin and J. Schwinger, Phys. Rev. 115, 1342 energies are calculated.
(1959). 42 A.F. Hebard et al., Nature (London) 350, 600 (1991).
14 L. Hedin, Phys. Rev. 139, A796 (1965). 43 O. Gunnarsson, Rev. Mod. Phys. 69, 575 (1997).
15 M.S.Hybertsen,S.G.Louie,Phys.Rev.B34,5390(1986). 44 V.P. Antropov, O. Gunnarsson, A.I. Liechtenstein, Phys.
16 R.W. Godby, M. Schlu¨ter, and L.J. Sham, Phys. Rev. B Rev. B 48, 7651 (1993).
37, 10159 (1988). 45 Coˆt´e,M.andGrossman,J.C.andCohen,M.L.andLouie,
17 G.Onida, L.Reining, A.Rubio, Rev.Mod.Phys.74, 601 S. G., Phys. Rev. Lett. 81, 697 (1998).
(2002). 46 N.Iwahara,T.Sato,K.Tanaka,andL.F.Chibotaru,Phys.
18 F.Aryasetiawan,O.Gunnarsson,Rep.Prog.Phys.61,237 Rev. B 82, 245409 (2010). For an interpretation of the
(1998). experimental data, see Note 21 of Ref. 9.
19 B. Farid, R. Daling, D. Lenstra, and W. van Haeringen, 47 I.D. Hands, J.L. Dunn, L. Janette, C.A. Bates, A. Colin,
Phys. Rev. B 38, 7530 (1988); B. Farid, in Electron Cor- M.J. Hope, J. Michael, S.R. Meech, D.L. Andrews, Phys.
relation in the Solid State, edited by N. H. March (World Rev. B 77, 115445 (2008). For an interpretation of the
Scientific,Singapore,1999),p.217,andreferencestherein. experimental data, see Note 22 of Ref. 9.
20 F.Bruneval,Ph.D.thesis,EcolePolytechnique,Palaiseau, 48 P.H.Hahn,W.G.Schmidt,andF.Bechstedt,Phys.Rev.
France, 2005. B72,245425(2005).Inthisstudy,theimportanceofcor-
21 F.Bruneval,N.Vast,L.Reining,Phys.Rev.B74,045102 rectingthestartingeigenvaluesintheconstructionofWis
(2006). emphasized for molecular systems, bringing some light in
22 M.Gatti,F.Bruneval,V.Olevano,L.Reining,Phys.Rev. thediscussionsconcerningtheupdateoftheeigenvaluesin
Lett. 99, 266402 (2007). theGreen’sfunctiononly(theGW scheme)orinbothG
0
23 F.Bruneval,N.Vast,L.Reining,M.Izquierdo,F.Sirotti, and W.
N. Barrett, Phys. Rev. Lett. 97, 267601 (2006). 49 L.Tiago,P.R.C.Kent,R.Q.Hood,andF.A.Reboredo,
24 F. Trani, J. Vidal, S. Botti, and M.A.L. Marques Phys. J. Chem. Phys. 129, 084311 (2008).
Rev. B 82, 085115 (2010). 50 Sahar Sharifzadeh, Ariel Biller, Leeor Kronik, and Jeffrey
25 J. Vidal, S. Botti, P. Olsson, J.-F. Guillemoles, and L. B. Neaton, Phys. Rev. B 85, 125307 (2012).
Reining, Phys. Rev. Lett. 104, 056401 (2010). 51 Conor Hogan, Maurizia Palummo, Johannes Gierschner,
26 M.J. Lucero, I. Aguilera, C. Diaconu, P. Palacios, P. and Angel Rubio, J. Chem. Phys. 138, 024312 (2013).
Wahno´n, and G. E. Scuseria, Phys. Rev. B 83, 205128 52 C. Faber, P. Boulanger, I. Duchemin, C. Attaccalite, and
(2011). X. Blase, J. Chem. Phys. 139, 194308 (2013).
27 W. Kang, M.S. Hybertsen, Phys. Rev. B 82, 195108 53 B. Baumeier, D. Andrienko, Y. Ma, and M. Rohlfing, J.
(2010). Chem. Theory Comput. 8, 997 (2012).
28 L.J. Sham and T.M. Rice, Phys. Rev. 144, 708 (1966). 54 P. Boulanger, D. Jacquemin, I. Duchemin, X. Blase, J.
29 W.HankeandL.J.Sham,Phys.Rev.Lett.43,387(1979). Chem. Theory Comput. 10, 1212 (2014). The effect of
30 G. Strinati, La Rivista del Nuovo Cimento 11 (12), 1-86 partial self-consistency on the eigenvalues, with various
(1988). starting-point wave functions, is discussed, in comparison
31 M. Rohlfing and S.G. Louie, Phys. Rev. Lett. 80, 3320 with high-level coupled-cluster data.
(1998). 55 A.Marini,C.Hogan,M.Gruu¨ning,andD.Varsano,Com-
32 L.X. Benedict, E.L. Shirley, and R.B. Bohn, Phys. Rev. put. Phys. Commun. 180, 1392 (2009).
Lett. 80, 4514 (1998). 56 X. Gonze et al., Comp. Mat. Science, 25, 478 (2002).
33 S.Albrecht,L.Reining,R.DelSole,G.Onida,Phys.Rev. 57 Notice that in Ref. 5 the interlayer distance between the
Lett. 80, 4510 (1998). graphene sheets is 15 a.u. and not 20 a.u. as it is wrongly
34 X.Blase,C.Attaccalite,V.Olevano,Phys.Rev.B83,536 reported in the paper.
115103 (2011). 58 P.B.AllenandV.Heine,J.Phys.C9,2305(1976);P.B.
35 C. Faber, C. Attaccalite, V. Olevano, E. Runge, X. Blase, Allen and M. Cardona, Phys. Rev. B 23, 1495 (1981).
Phys. Rev. B 83, 115123 (2011). 59 S. Ponc´e, G. Antonius, P. Boulanger, et al; Comput.
36 X.BlaseandC.Attaccalite,Appl.Phys.Lett.99, 171909 Mater. Sci. 83, 341 (2014).
(2011). 60 For sake of indication, at the LDA level, the two lowest
37 J.M.Soler,E.Artacho,J.D.Gale,A.Garc´ıa,J.Junquera, X and X vibrational modes lead to a (negative) -0.14
1 2
P.Ordejo´n,D.Sa´nchez-Portal,J.Phys.: Condens.Matter eV coupling energy, one order of magnitude smaller than
14, 2745 (2002). the X mode coupling strength.
4
38 Thistriple-zetapluspolarization(TZP)basisconsistsof3 61 P. Trevisanutto, C. Giorgetti, L. Reining, M. Ladisa and
shells of s and p orbitals, complemented by one polariza- V. Olevano, Phys. Rev. Lett. 101, 226405, (2008)
tion d shell, that is 19 orbitals per carbon atom, for the 62 We provide for information in Table II the results asso-
10
ciated with the (single-shot) COHSEX approach, starting Marques, S. Botti, J. Chem. Theory Comput. 10, 3934
fromDFT-LDAeigenstates.Clearly,theEPCcouplingen- (2014).
ergiesapproachfrombelowtheevGW valuewithanerror 67 TheexperimentalbandgapistakenfromtheNISTchem-
notlargerthanthatinducedbytheevCOHSEXtechnique. istry webbook at [http://webbook.nist.gov/chemistry/].
ThisstemsfromacancellationoferrorsbetweentheCOH- Thebandgapiscalculatedasthedifferenceoftheexperi-
SEXapproximationandthelackofself-consistencywhich mental ionization energy and electronic affinity. Weadopt
leads to overscreening when building W directly from the the most recent experimental data.
DFT-LDA spectrum. 68 DirectcomparisonwithFigure2ofRef.12isnotpossible,
63 S. Ismail-Beigi and S. G. Louie, Phys. Rev. Lett. 90, sincethefigurepresentsthetotalelectron-phononcoupling
076401 (2003). atwavevectorq=X,i.e.allphononmodesareconsidered
64 Thomas Ko¨rzdo¨rfer and Noa Marom, Phys. Rev. B 86, andnotjustX asinthisstudy.Thecomparisonwasmade
4
041110(R) (2012). directly with the data used in Ref. 12 via personal corre-
65 F. Bruneval, M.A.L. Marques, J. Chem. Theory Comput. spondence.
9, 324 (2013).
66 S. Korbel, P. Boulanger, I. Duchemin, X. Blase, M.A.L.
