Large scale GW calculations

5.7 MB·English

by  Marco Govoni

Large scale GW calculations
Marco Govoni ,†,‡ and Giulia Galli ,†,‡
∗ ∗
5
1
0
InstituteforMolecularEngineering,TheUniversityofChicago,ChicagoIL,USA,andMaterials
2

n
ScienceDivision,ArgonneNationalLaboratory,ArgonneIL,USA
a
J

3
1 E-mail: [email protected];[email protected]

]
i
c
s
-
l Abstract
r
t
m
We present GW calculations of molecules, ordered and disordered solids and interfaces,
.
t
a
m which employ an efficient contour deformation technique for frequency integration, and do
-
d not require the explicit evaluation of virtual electronic states, nor the inversion of dielectric
n
o matrices. Wealsopresentaparallelimplementationofthealgorithmwhichtakesadvantageof
c
[
separable expressions of both the single particle Green’s function and the screened Coulomb

1
v interaction. The method can be used starting from density functional theory calculations per-
1
4
formed with semi-local or hybrid functionals. We applied the newly developed technique to
1
3
GW calculations of systems of unprecedented size, including water/semiconductor interfaces
0
.
1
withthousandsofelectrons.
0
5
1
:
v
i 1 Introduction
X
r
a
Theaccuratedescriptionoftheexcitedstatepropertiesofelectronsplaysanimportantroleinmany
fieldsofchemistry,physics,andmaterialsscience1. Forexample,theinterpretationandprediction
Towhomcorrespondenceshouldbeaddressed
∗
†InstituteforMolecularEngineering,TheUniversityofChicago,ChicagoIL,USA
‡MaterialsScienceDivision,ArgonneNationalLaboratory,ArgonneIL,USA
1
of photoemission and opto-electronic spectra of molecules and solids rely on the ability to com-
putetransitionsbetweenoccupiedandvirtualelectronicstatesfromfirstprinciples,aswellastheir
lifetimes2.
In particular, in the growing field of materials for energy conversion processes – including solar
energy conversion by the photovoltaic effect and solar to fuel generation by water photocatalysis
–ithasbecomekeytodeveloppredictivetoolstoinvestigatetheexcitedstatepropertiesofnanos-
tructures and nanocomposites and of complex interfaces3–5. The latter include solid/solid and
solid/liquidinterfaces,e.g. betweenasemiconductororinsulatorandwateroranelectrolyte6–10.
In the last three decades, Density Functional Theory (DFT) has been widely used to compute
structural and electronic properties of solids and molecules11–15. Although successful in describ-
ing ground state and thermodynamic properties, and in many ab initio molecular dynamics stud-
ies16,17,DFTwithbothsemi-localandhybridfunctionalshasfailedtoaccuratelydescribeexcited
state properties of several materials18. However, hybrid functionals have brought great improve-
ment for properties computed with semi-local ones, e.g. for defects in semiconductor and ox-
ides19–22. In particular hybrid functionals with admixing parameters computed self-consistently
have shown good performance in reproducing experimental band gaps and dielectric constants of
broad classes of systems23. In the case of the electronic properties of surfaces, interfaces (and
hence nanostructures), the use of hybrid functionals has in many instances not been satisfactory.
Indeed calculations with hybrid functionals yield results for electronic levels that often depend on
the mixing parameter chosen for the Hartree-Fock exchange; such parameter is system dependent
andthereisnoknownfunctionalyieldingsatisfactoryresultsfortheelectronicpropertiesofinter-
faces composed of materials with substantially different dielectric properties, as different as those
of, e.g. water (ε = 1.78)24 and Si (ε = 11.9)25 or water and transition metal oxides of interest
∞ ∞
forphotoelectrodes(ε =5-7)26.
∞
The use of many body perturbation (MBPT) starting from DFT single particle states has proven
accurateforseveralclassesofsystems27–36 anditappearstobeapromisingframeworktodescribe
complex nanocomposites and heterogeneous interfaces. MBPT for the calculations of photoemis-
2
sion spectra in the GW approximation37, and of optical spectra by solving approximate forms of
the Bethe Salpeter Equation (BSE)38 is in principle of general applicability; however its use has
beengreatlylimitedbycomputationaldifficultiesinsolvingtheDyson’sequationandtheBSEfor
realisticsystems.
Recently we proposed a method to compute quasi particle energies within the G W approxima-
0 0
tion(i.e. thenon-selfconsistentGW approximation)thatdoesnotrequiretheexplicitcalculationof
virtual electronic states, nor the inversion of large dielectric matrices39,40. In addition the method
does not use plasmon pole models but instead frequency integrations are explicitly performed and
there is one single parameter that controls the accuracy of the computed energies, i.e. the number
of eigenvectors and eigenvalues used in the spectral decomposition of the dielectric matrix at zero
frequency. The method was successfully used to compute the electronic properties of water41 and
aqueoussolutions42andofheterogeneoussolids5,includingcrystallineandamorphoussamples40.
However the original method contained some numerical approximations in the calculations of the
headandwingsofthepolarizabilitymatrix;mostimportantlythecorrelationself-energywascom-
puted on the imaginary axis and obtained in real space by analytic continuation. Finally, although
exhibitingexcellent scalability, themethod wasnotyetapplied tosystemswith thousandsofelec-
trons,e.g. torealisticinterfaces,duetothelackofparallelizationinitsoriginalimplementation.
In this paper we solved all of the problems listed above, by (i) eliminating numerical approxima-
tions in the calculation of the polarizability; (ii) avoiding the use of an analytic continuation and
using efficient contour deformation techniques; (iii) providing a parallel implementation of the
algorithm based on separable forms of both the single particle Green function and the screened
Coulombinteraction. ThemethodpresentedherecanbeusedstartingfromDFTorbitalsandener-
giesobtainedbothwithsemi-localandhybridfunctionals. Weappliedourtechniquetothecalcula-
tion of the electronic properties of systems of unprecedented size, including water/semiconductor
interfaceswithmorethanonethousandelectrons. Thesecalculationsallowonetoaccuratelystudy
states at heterogeneous interfaces and to define an electronic thickness of solid/liquid interfaces
usingMBPT.
3
Therestofthepaperisorganizedasfollows. Sec.Section2describestheG W methodologythat
0 0
we implemented in a computational package called West. Sec. Section 3 presents results for the
ionization potentials of closed and open shell molecules and for the electronic structure of crys-
talline, amorphous and liquid systems, aimed at verifying and validating the algorithm and code
West against previous calculations and measurements. Sec. Section 4 presents the study of the
electronic properties of finite and extended large systems, i.e. nanocrystals and solid/liquid inter-
faces,ofinteresttophotovoltaicandphotocatalysisapplications,respectively. Ourconclusionsare
reportedinSec.Section5.
4
2 Method
WithinDFT,then-thsingleparticlestateψ andenergyε ofasystemofinteractingelectrons
nkσ nkσ
isobtainedbysolvingtheKohn-Sham(KS)equation11–15:
Hˆσ ψ =ε ψ (1)
KS nkσ nkσ nkσ
| (cid:105) | (cid:105)
where Hˆσ = Tˆ +Vˆ +Vˆ +Vˆσ is the KS Hamiltonian, Tˆ is the kinetic energy operator, and
KS ion H xc
Vˆ ,Vˆ andVˆσ are the ionic, Hartree, and exchange-correlation potential operators, respectively.
ion H xc
The indexes k and σ label a wavevector within the first Brillouin zone (BZ) and spin polarization,
respectively. Hereweconsidercollinearspin,i.e. decoupledupanddownspins.
QP
InafashionsimilartoEq.(Eq.(1))onemayobtainquasiparticle(QP)statesψ andQPenergies
nkσ
QP
E bysolvingtheequation:
nkσ
Hˆσ ψQP =EQP ψQP (2)
QP nkσ nkσ nkσ
(cid:12) (cid:69) (cid:12) (cid:69)
(cid:12) (cid:12)
where the QP Hamiltonian Hˆσ is form(cid:12)ally obtained by(cid:12)replacing, in Eq. (Eq.(1)), the exchange-
QP
correlation potential operator with the electron self-energy operator Σσ = iGσWΓ; Gσ is the in-
teractingone-particleGreen’sfunction,W isthescreenedCoulombinteractionandΓisthevertex
operator28,43. All quantities entering the definition of the self-energy are interdependent and can
be obtained self-consistently adopting the scheme suggested by L. Hedin44–46. In the GW ap-
proximation, Γ is set equal to the identity, which yields the following expression for the electron
self-energy47:
+∞
dω
Σσ(r,r;ω)=i (cid:48)Gσ(r,r;ω+ω )W (r;r;ω ), (3)
(cid:48) (cid:48) (cid:48) RPA (cid:48) (cid:48)
2π
(cid:90)
∞
−
where W is the screened Coulomb interaction obtained in the random phase approximation
RPA
(RPA). Due to the non-locality and the frequency dependence of the electron self-energy, a self-
consistent solution of Eq. (Eq.(2)) is computationally very demanding also for relatively small
5
QP
systems,containingtensofelectrons,andusuallyoneevaluatesQPenergiesE perturbatively:
nkσ
EQP = ε + ψ Hˆσ Hˆσ ψ (4)
nkσ nkσ (cid:104) nkσ| QP− KS | nkσ(cid:105)
= ε + ψ Σ(cid:0)ˆσ(EQP ) ψ(cid:1) ψ Vˆσ ψ . (5)
nkσ (cid:104) nkσ| nkσ | nkσ(cid:105)−(cid:104) nkσ| xc| nkσ(cid:105)
QP
We note that E enters both the left and right hand side of Eq. (Eq.(3)), whose solution is
nkσ
usually obtained recursively, e.g. with root-finding algorithms such as the secant method. The use
ofEq.(Eq.(3))toevaluateQPenergiesfromKSstatesandofthecorrespondingKSwavefunctions
isknownastheG W approximation.
0 0
WithinG W ,usingtheLehmann’srepresentation,theGreen’sfunctionis:
0 0
dk ψ (r)ψ (r)
Gσ (r,r;ω)= ∑ nkσ n∗kσ (cid:48) (6)
KS (cid:48) − (2π)3ε ω iηsign(ε ε )
nB(cid:90)Z nkσ − − nkσ − F
where η is a small positive quantity and ε is the Fermi energy. In Eq. (Eq.(6)) we used the
F
subscript KS to indicate that the Green’s function is evaluated using the KS orbitals obtained by
solvingEq.(Eq.(1)).
InRef.[39,40]analgorithmwasintroducedtocomputetheself-energymatrixelementsofEq.(Eq.(3))
withouttheneedtoevaluateexplicitlyempty(virtual)electronicstates,byusingatechniquecalled
projective eigendecomposition of the dielectric screening (PDEP). A diagram of the method is re-
portedinFig.Figure1. AfterKSsingleparticleorbitalsandenergiesareobtainedusingsemilocal
or hybrid functionals, the screened Coulomb interaction is computed using a basis set built from
the eigenpotentials of the static dielectric matrix at zero frequency. In this way W entering
RPA
Eq.(Eq.(3))isexpressedinaseparableform,similartothatofGσ inEq.(Eq.(6)). Inthefollow-
KS
ingsectionswedescribeindetailallthestepsoutlinedinFig.Figure1. TheseparableformofW
RPA
isgiveninSec.Section2.1. Calculationofthepolarizabilityandthespectraldecompositionofthe
dielectricmatrixaredescribedinSec.Section2.2andSection2.3,respectively. Matrixelementsof
G andW are then obtained without the explicit use of empty electronic states and simultane-
KS RPA
ously at several frequencies by using a deflated Lanczos technique, described in Sec. Section 2.4.
6
Finally the frequency integration is carried out by introducing a contour deformation method, as
described in Sec. Section 2.5. The use of the analytic continuation used in the original method of
Ref.[39,40]isthusavoided.
DDFFTT PDEP Lanczos W
Contour
Lanczos G
Deformation
⌃(!)= d!0G(!+!0)W(!0)
Z
QP-energies
Figure 1: (Color online) Schematic representation of the steps involved in the calculations of
quasiparticle (QP) energies, within the G W approximation, using the method proposed in this
0 0
work. The KS energies (ε) and occupied orbitals (ψ) computed at the DFT level are input to the
i i
PDEP algorithm, which is used to iteratively diagonalize the static dielectric matrix (ε 1) at zero
−
frequency. Thesetofeigenvectors φ constitutesthebasissetusedtocomputebothGandW at
n
{ }
finitefrequencieswiththeLanczosalgorithm. ThefrequencyintegrationofEq.(Eq.(3))iscarried
out using the contour deformation technique. The frequency dependent matrix elements of the
electron self-energy are thus obtained and introduced in Eq. (Eq.(3)) to compute the QP energies
QP
E .
i
2.1 Separable form of the screened Coulomb interaction
In order to solve equation Eq. (Eq.(3)) and obtain QP energies, one needs to compute the ma-
trix elements of the electron self-energy between KS states, which in the G W approximation is
0 0
given by Eq. (Eq.(3)). The Green’s function may be expressed in a fully separable form using its
Lehmann’srepresentation,Eq.(Eq.(6)). HoweverW isnottriviallyseparable;itisgivenasthe
RPA
sumoftwoterms:
W (r,r;ω)=v(r,r)+W (r,r;ω) (7)
RPA (cid:48) (cid:48) p (cid:48)
7
where v(r,r)= e2 is the bare Coulomb interaction48 and W is a nonlocal and frequency de-
(cid:48) r r p
| − (cid:48)|
pendent function. Using Eq. (Eq.(7)), we write the self-energy as the sum of two contributions
Σσ =Σσ +Σσ,whereonlythelatterdependsonthefrequency:
X C
+∞
dω
Σσ(r,r) = i (cid:48)Gσ (r,r;ω+ω )v(r;r) (8)
X (cid:48) 2π KS (cid:48) (cid:48) (cid:48)
(cid:90)
∞
−
Nσ
occ dk
= ∑ ψ (r)v(r,r)ψ (r) (9)
− (2π)3 nkσ (cid:48) n∗kσ (cid:48)
n=1(cid:90)
BZ
Nσ is the number of occupied states with spin σ; Σσ is usually called exchange self-energy
occ X
becauseitisformallyequivalenttotheFockexactexchangeoperator49;
+∞
dω
Σσ(r,r;ω)=i (cid:48)Gσ (r,r;ω+ω )W (r,r;ω ) (10)
C (cid:48) 2π KS (cid:48) (cid:48) p (cid:48) (cid:48)
(cid:90)
∞
−
Σσ is referred to as correlation self-energy. Using Eq.s (Eq.(7))-(Eq.(10)) the QP Hamiltonian of
C
Eq.(Eq.(2))maybeexpressedas:
Hˆ (ω)=Tˆ +Vˆ +Vˆσ +Σˆσ(ω) (11)
QP ion HF C
where Vˆσ is the Hartree-Fock potential operator. The ionic potential Vˆ is treated within the
HF ion
pseudopotentialapproach50.
In this work we expressW in a separable form by adopting the projective dielectric eigendecom-
p
position (PDEP) technique, proposed in Ref. [51-52], and we use a plane wave basis set to express
thesingleparticlewavefunctionsandchargedensity,withinperiodicboundaryconditions:
ψ (r)=eikru (r)=∑c (G)ei(k+G)r (12)
nkσ · nkσ nkσ ·
G
where G is a reciprocal lattice vector, c (G) = 1 dru (r)e iGr and Ω is the unit cell
nkσ Ω Ω nkσ − ·
volume. In Eq. (Eq.(12)) all reciprocal lattice vectors (cid:82)such that 1 k+G 2 <E are included
2| | cutwfc
8
inthesummation. UsingaplanewavedescriptionalsoforW wehave
p
dq
W (r,r;ω)= ∑ ei(q+G)r v +Wp (q;ω) e i(q+G(cid:48))r(cid:48) (13)
RPA (cid:48) (2π)3 · GG(cid:48) GG(cid:48) − ·
(cid:90) GG
BZ (cid:48)
(cid:2) (cid:3)
wherev = 4πe2 δ (δ istheKroneckerdelta)and
GG(cid:48) q+G2 GG(cid:48)
| |
χ¯ (q;ω)
p GG
W (q;ω)= (cid:48) . (14)
GG(cid:48) q+G q+G
(cid:48)
| || |
In Eq. (Eq.(14)) we have introduced the symmetrized reducible polarizability χ¯, related to the
symmetrizedinversedielectricmatrixε¯ 1 bytherelation:
−
ε¯ 1 (q;ω)=δ +χ¯ (q;ω). (15)
G−G(cid:48) GG(cid:48) GG(cid:48)
Thesymmetrizedform χ¯ ofthepolarizability χ is
√4πe2 √4πe2
χ¯ = χ . (16)
GG GG
(cid:48) q+G (cid:48) q+G
(cid:48)
| | | |
The reducible polarizability χ is related to the irreducible polarizability χ by the Dyson’s equa-
0
tion,whichwithintheRPAreads:
χ =χ0 + ∑ χ0 v χ (17)
GG(cid:48) GG(cid:48) GG1 G1G2 G2G(cid:48)
G1,G2
orintermsofsymmetrizedpolarizabilites:
χ¯ =(1 χ¯0) 1χ¯0. (18)
−
−
Within a plane wave representation each quantity in Eq. (Eq.(18)) is a matrix of dimension N2 ,
pw
and in principle χ¯ can be obtained from χ¯0 via simple linear algebra operations. In practice,
the matrices of Eq. (Eq.(18)) may contain millions of rows and columns for realistic systems;
9
for example for a silicon nanocrystal with 465 atoms, placed in a cubic box of edge 90bohr, 1.5
million plane waves are needed in the expansion of Eq. (Eq.(12)) with E =25Ry. It is thus
cutwfc
important to find alternative representations of χ¯ and reduce the number of elements to compute.
Onecouldthinkofastraightforwardspectraldecomposition:
N
pdep
χ¯ (q;ω)= ∑ φ (q+G;ω)λ (q;ω)φ (q+G;ω) (19)
GG i i ∗j
(cid:48)
i=1
where φ and λ are the eigenvectors and eignvalues of χ¯, respectively. Unfortunately this strategy
i i
is still too demanding from a computational standpoint, as it implies finding eigenvectors and
eigenvaluesatmultiplefrequencies.
Acomputationallymoretractablerepresentationmaybeobtainedusingthespectraldecomposition
of χ¯ at ω = 0. As apparent from Eq. (Eq.(18)), eigenvectors of χ¯ are also eigenvectors of χ¯0;
0
the latter is easier to iteratively diagonalize than χ¯, and the frequency dependency may be dealt
with iterative techniques, starting from the solution at ω = 0, as discussed in Sec. Section 2.4.
Therefore we proceed by solving the secular equation for χ¯0 only at ω =0, generating what we
call the PDEP basis set φ : i=1,N , which is used throughout this work to represent the
i pdep
| (cid:105)
polarizability χ¯: (cid:8) (cid:9)
N
pdep
χ¯ (q;ω)= ∑ φ (q+G)Λ (q;ω)φ (q+G); (20)
GG i ij ∗j
(cid:48)
i=1,
j=1
hereΛ (q;ω)isamatrixofdimensionN2 . IngeneralN N ,51,52 leadingtosubstantial
ij pdep pdep pw
(cid:28)
computational savings.53 The N functions φ may be computed by solving the Sternheimer
pdep i
equation54,withoutexplicitlyevaluatingempty(virtual)electronicstates. Inaddition,N turns
pdep
out to be the only parameter that controls the accuracy of the expansion in Eq. (Eq.(20)). The
detailsofthederivationofthePDEPbasissetaregiveninSec.Section2.3. Wenotethatalternative
basis sets, based on the concepts related to maximally localized Wannier functions, have been
proposedintheliteraturetoreducethedimensionalityofthepolarizabilitymatrices55.
10