Gaussian-based coupled-cluster theory for the ground state and band structure of solids

0.34 MB·

by  James McClain

Gaussian-based coupled-cluster theory for the ground state and band
structure of solids
JamesMcClain,1 QimingSun,2 GarnetKin-LicChan,2,a) andTimothyC.Berkelbach3,b)
1)DepartmentofChemistry,PrincetonUniversity,Princeton,NewJersey08544,USA
2)DivisionofChemistryandChemicalEngineering,CaliforniaInstituteofTechnology,Pasadena,California91125,
USA
3)Department of Chemistry and James Franck Institute, University of Chicago, Chicago, Illinois 60637,
USA
WepresenttheresultsofGaussian-basedground-stateandexcited-stateequation-of-motioncoupled-clustertheorywith
singleanddoubleexcitationsforthree-dimensionalsolids. Wefocusondiamondandsilicon,whichareparadigmatic
7 covalentsemiconductors. Inaddition toground-state properties(the latticeconstant, bulkmodulus, andcohesive en-
1 ergy),wecomputethequasiparticlebandstructureandbandgap. WesampletheBrillouinzonewithupto64k-points
0 usingnorm-conservingpseudopotentialsandpolarizeddouble-andtriple-zetabasissets,leadingtocanonicalcoupled-
2 clustercalculationswithasmanyas256electronsin2,176orbitals.

n
a
J I. INTRODUCTION basedapproachesforsolids,namelyincrementalmethods26–29
7 andembeddingtechniques.30–32
1
The electronic ground-state and excited-state structure of
three-dimensional solids and low-dimensional nanomaterials Recently, we performed the first excited-state CC treat-
]
i isanongoingchallengeincomputationalmaterialsscience.A ment of a three-dimensional condensed-phase system, at the
c
theoreticalprogrambasedondensityfunctionaltheory(DFT), equation-of-motion (EOM)-CCSD and EOM-CCSDT level,
s
- frequentlycombinedwithtime-dependentmany-bodypertur- focusingontheuniformelectrongasatthedensityofmetallic
rl bation theory, has become a reliable and accurate standard sodium.33 Forsmallunitcellsinafinitesingle-particlebasis,
mt approach. In particular, ground-state properties of weakly theEOM-CCSDTone-particlespectralfunctionwasinnearly
correlatedmaterialsarewell-describedbyDFT,theprincipal perfect agreement with numerically exact (time-dependent
.
t chargedexcitedstatesleadingtobandstructurecanbetreated densitymatrixrenormalizationgroup)results,andmoreaccu-
a
m withhybrid1–3orsemilocal4functionals,aswellasthemany- ratethananytestedflavoroftheGWapproximation,whilethe
bodyGW approximation,5,6 andneutralexcitedstatescanbe EOM-CCSDone-particlespectralfunctionquantitativelyim-
d- described with time-dependent DFT7,8 or the Bethe-Salpeter proved on the standardGW approximation. Importantly, un-
n equation.9 However,errorsinthisapproachcanbedifficultto liketheGW approximation(inpractice,G0W0),theEOM-CC
o assess, due to (for example) the underlying DFT exchange- spectra are invariant to the starting choice of single-particle
c correlationfunctionalorself-consistencyintheGW approxi- energies. At larger system sizes approaching the thermody-
[
mation. Evenintheformalmany-bodyframeworks, system- namiclimit,itwasnotpossibletocarryoutfullEOM-CCSDT

1 aticimprovements(vertexcorrections)areonlyrecentlybeing calculations, but it was still possible to compute the EOM-
v pursued,withmixedsuccess.10–13 CCSDspectrum. SimilarlytotheGW-plus-cumulantapprox-
2 In contrast, wavefunction-based approaches originating in imation,34,35theEOM-CCSDspectrumremovedthefictitious
3 thequantumchemistrycommunityarenowstartingtobeap- plasmaronpoleoftheGWspectrum,36,37demonstratingthatit
8
plied to the properties of condensed-phase materials. Atom- correctlytreatsthecouplingofchargedexcitationstobosonic
4
0 isticone-andtwo-dimensionalsystemshavebeenstudiedby plasmonexcitationsinmetals.
. configuration interaction with single excitations (CIS), time-
1
dependent Hartree-Fock (HF), and time-dependent DFT,14 Here,wecontinuethisprogramandpresentthefirstresults
0
7 as well as second-order Møller-Plesset perturbation the- of CC excited states for atomistic solids. The layout of this
1 ory (MP2)15–18 and coupled-cluster (CC) theory with single paper is as follows. In Sec. II, we describe the Gaussian ba-
: anddoubleexcitations(CCSD).15,17,18 Forthree-dimensional sis sets, integral evaluation, and HF calculations for periodic
v
solids, MP2 has been implemented with local correlation in systems,includingissuesassociatedwithdivergenttermsand
i
X the CRYSCOR program,19,20 in the plane-wave based VASP finite-size effects. In Sec. III, we describe the ground-state
r package with full Brillouin zone sampling,21–23 and in the CCSDandexcited-stateEOM-CCSDformalismsforperiodic
a mixed Gaussian and plane-wave CP2K package at the Γ- systems,aswellastheconnectiontoquasiparticlebandstruc-
point.24 Most recently, ground-state CCSD, CCSD with per- ture. In Sec. IV we discuss the convergence and finite-size
turbative triple excitations (CCSD(T)), and full CI quantum effects of our HF and correlated calculations. We present
Monte Carlo for periodic solids have been reported.23,25 We ground-state HF, MP2, and CCSD results in Sec. V for di-
alsomentiontheparalleldevelopmentofnon-periodiccluster- amond and silicon. This section also includes the techni-
cal details of our calculations, a discussion of convergence,
and comparison to existing MP2 and CCSD calculations on
solids.22,25 In Sec. VI, we present EOM-CCSD band struc-
a)Electronicmail:[email protected] tures and a convergence study of the indirect band gap. We
b)Electronicmail:[email protected] concludeinSec.VII.
2
II. GAUSSIAN-BASEDPERIODICINTEGRALSAND contributionfromtheionpseudopotentialsvL(r)isgivenby,
HARTREE-FOCKTHEORY (cid:88)(cid:88)
vL(r)= eiG·(r−RI)vLI(G), (6)
We use an underlying single-particle basis of crys- G(cid:44)0 I
talline Gaussian-based atomic orbitals (AOs). These are where I indexes the ions at position R in the unit cell. The
I
translational-symmetry-adaptedlinearcombinationsofGaus- local part of the GTH pseudopotential is finite at r = 0 but
sianAOsoftheform decays as Z /r leading to a divergent G = 0 component,
I
(cid:88) which is tr−eated separately (see below). The nonlocal (NL)
φµk(r)= eik·Tφ˜µ(r−T)≡eik·ruµk(r) (1) pseudopotentialpartisevaluatedonthereciprocal-spacegrid
T
(cid:88)
where T is a latice translation vector and k is a crystal mo- VµNνL(k)=Ω φ∗µk(G)vNL(k+G,k+G(cid:48))φνk(G(cid:48)) (7)
mentum vector in the first Brillouin zone. In this work, we G,G
(cid:48)
sample k from a uniform (but not necessarily Monkhorst-
with
Pack)grid. ThesecondequalityaboveexpressesBloch’sthe-
(cid:88)
orem, wheretheBlochfunctionuµk(r)isfullyperiodicwith vNL(k+G,k+G(cid:48))= e−i(G−G(cid:48))·RIvNIL(k+G,k+G(cid:48)) (8)
respect to all lattice translations. Therefore, we can exactly I
expandthecrystallineAOsinasetofauxiliaryplane-waves
and
(cid:88)
φµk(r)= G(cid:90)φµk(G)ei(k+G)·r, (2a) vNIL(k+G,k+G(cid:48))= Ω1 (cid:90) dr(cid:90) dr(cid:48)e−i(k+G)·r (9)
φµk(G)= Ω1 drφµk(r)e−i(k+G)·r, (2b) ×vNIL(r,r(cid:48))ei(k+G(cid:48))·r(cid:48).
Ω
In the above, the ion-specific local and nonlocal parts of
where G is a reciprocal lattice vector and Ω is the unit cell
the pseudopotential have analytic definitions of the GTH
volume.
form.40,41 In particular the nonlocal part is separable and ex-
Giventhisbasischoice, thereareseveraldifferentwaysto
pressibleasasumofproductsoffunctionsofk+Gandfunc-
evaluate the corresponding overlap, kinetic energy, electron-
tionsofk+G.
(cid:48)
nuclearattraction,andCoulombintegralsarisinginquantum
The Hartree and exchange matrices are evaluated in real
chemistry. We will present a careful comparison of different
space,
choices in a later publication, and here only detail the grid-
(cid:90)
based scheme that we have used for the results in this work,
whichcloselyresemblestheGaussianandplane-wavescheme Jµν(k)= drφ∗µk(r)vH(r)φνk(r), (10)
Ω
usedintheCP2Kpackage.38,39 Weuseadualreal-spaceand (cid:90) (cid:90) ρ(r,r )
reciprocal-spacerepresentationofthecrystallineAOs. Inreal Kµν(k)= dr dr(cid:48)φ∗µk(r) r r(cid:48) φνk(r(cid:48)), (11)
space, we represent φµk(r) on a uniform real-space grid in Ω | − (cid:48)|
the unit cell. In reciprocal space, we obtain φ (G) by a intermsoftheHartreepotential
µk
fastFouriertransform(FFT)leadingtoarepresentationona
4π (cid:88) ρ(G)
uniformgridofpointscorrespondingtothereciprocallattice vH(r)= Ω G2 eiG·r (12)
vectors of our unit cell. The real-space grid density roughly G(cid:44)0
correspondstoakineticenergycutoffintheplane-waverep-
andthedensitymatrixρ(r,r ). Ingeneral,thedensitymatrix
resentation. Toensurethatrelativelylowcutoffscanbeused, (cid:48)
and density ρ(r) = ρ(r,r) can be obtained from a Brillouin
wereplacethecoreelectronswithseparablenorm-conserving
zonesampling
GTH (HGH) pseudopotentials,40,41 which removes the sharp
nucleardensities. WeusetheGaussianbasissetsthatarede- 1 (cid:88)(cid:88)
signed for use with these pseudopotentials in solid-state cal- ρ(r,r(cid:48))= N Pλσ(k)φλk(r)φ∗σk(r(cid:48)), (13)
culations.38 k k λσ
The one-electron overlap, kinetic energy, and local (L) where N is the number of k-points sampled in the Brillouin
k
pseudopotentialintegralsareevaluatedbynumericalintegra- zone. Assuming a closed-shell reference of molecular or-
(cid:80)
tiononthereal-spacegridaccordingto bitals, ψ = C (k)φ , the mean-field density matrix
pk µ µp µk
(cid:90) isgivenby
Sµν(k)= drφ∗µk(r)φνk(r), (3) (cid:88)occ
Ω
1(cid:90) Pλσ(k)=2 Cλi(k)Cσ∗i(k). (14)
Tµν(k)=−2 drφ∗µk(r)∇2rφνk(r), (4) i
Ω
(cid:90) Thisleadstothereal-spaceformoftheexchangematrix
where we nVoµLtνe(kt)ha=t tΩhedsreφ∗µinkt(erg)rvaLl(sr)aφnνdk(ro)t;her quantit(i5es) Kµν(k)= N1k (cid:88)k (cid:88)λσ Pλσ(k(cid:48))(cid:90) drφ∗µk(r)vXνk,σk(cid:48)(r)φλk(cid:48)(r)
(cid:48)
throughoutthepaperaredefinedperunitcell. Thetotal,local (15)
3
where where the primed summation neglects the self-interaction
vXσk(cid:48),νk(r)=(cid:90) dr(cid:48)φ∗σk(cid:48)|(rr−(cid:48))φrν(cid:48)k|(r(cid:48)) ≡(cid:90) dr(cid:48)ρσ|rk(cid:48)−,νkr(r(cid:48)|(cid:48)) tteiorFmnossriswacitshhuobIsse=enqtJuoewgnehttehcneorrTrweli=atht0ioηantnotdrfeathacteimlirtaeanntegt,erwaopefidlnaectteoidcnevteosrugdmeenmficnaee-.
= 4π(cid:88)(cid:48) ρσk(cid:48),νk(G) ei(k−k(cid:48)+G)·r, themolecularorbitalintegrals. Theone-electronintegralscan
Ω k k +G2 be obtained straightforwardly by changing basis with C(k).
G | − (cid:48) |
Two-electronintegralsaredefined(again,perunitcell)as
(16)
andtheG =0termisexcludedwhenk=k(cid:48).Liketheindivid- (cid:104)pkp(cid:90),rkr|qk(cid:90)q,sks(cid:105)=(pkp,qkq|rkr,sks)
ualcrystallineAOs, theAOpairdensitieshaveaplane-wave
resolution, = Ωdr1 dr2ψ∗pkp(r1)ψqkq(r1)v(r12)ψ∗rkr(r2)ψsks(r2).
ρσk(cid:48),νk(r)=φ(cid:88)∗σk(cid:48)(r)φνk(r) (23)
= ρσk,νk(G)ei(k−k(cid:48)+G)·r, (17a) Inthepresenceoftranslationalinvariance,thesetwo-electron
(cid:48)
G integralsmustconservecrystalmomentum,i.e.kp+kr kq
(cid:90) − −
1 k = G, where G is a reciprocal lattice vector. We evalu-
ρσk(cid:48),νk(G)= Ω ρσk(cid:48),νk(r)e−i(k−k(cid:48)+G)·r. (17b) atse these integrals by representing the molecular orbital pair
Ω
densityintermsofplane-waves,leadingto
As shown above, the local part of the pseudopotential and
Hartree potential are separately divergent atG = 0, however (cid:88)(cid:48) (cid:104)
theirsumisnot. Foreachatom,thesumisgivenby (pkp,qkq|rkr,sks)=Ω2 ρpkp,qkq(G)
(cid:90) (cid:32) (cid:33) G (24)
Z e2 (cid:105)
αI = dr vLI(r)+ Ir , (18) ×v(kq−kp+G)ρrkr,sks(Gprqs−G) .
where αI is a finite, ion-specific parameter of the pseudopo- whereGprqs =kp+kr−kq−ksisareciprocallatticevector,
tential. Thisleadstoaadditionalmatrixelement v(G) = 4π/ΩG2,andtheG = 0singularityisremovedwhen
VµLν+J(k)= SµΩν(k)(cid:88)αI, (19) kpO=nckeq.theseintegralsareobtained,wecandefinethesecond-
I quantized Hamiltonian suitable for a correlated electronic
resultingsimplyinauniformshiftoftheorbitalenergiesanda structure treatment. For closed-shell single-reference corre-
constantinthetotalenergy. Notethatintheabsenceofpseu- lationtechniques,itisconvenienttoexpresstheHamiltonian
dopotentials,thedivergentcontributionsfromtheHartreeand innormal-orderedform
electron-nuclearinteractioncancelexactlybychargeneutral- H = ΦHΦ +(cid:88)(cid:88)F (k)(cid:110)Epk(cid:111)
ity. (cid:104) | | (cid:105) pq qk
pq k
Oncetheaboveintegralsaredefined, thenitisstraightfor-
wardtocarryoutaHFcalculationusingtheseintegrals. The + 1(cid:88) (cid:88)(cid:48) (pk qk rk sk )(cid:26)EpkpErkr(cid:27) (25)
onlydifferencecomparedtoastandardmolecularHFcalcula- 2 p q| r s qkq sks
pqrsk kkk
tionisthattheintegralsandorbitalsarecomplex. Themolec- p q r s
euqlauratoiorbnitals at each k-point ψpk are obtained from the HF where Eqpkkqp = (cid:80)σa†pkpσaqkqσ is a spin-summed excitation
operator, ... denotes normal ordering with respect to the
F(k)C(k)=ε(k)S(k)C(k) (20) reference{Φ (}the Fermi vacuum), and as usual the primed
summation enforces crystal momentum conservation. Natu-
where
rally,inthecaseofcanonicalHF,theFockmatrixisdiagonal,
F(k)=T(k)+VPP(k)+J(k) 1K(k)+VL+J(k) (21) F (k)=ε (k)δ .
− 2 pq p pq
and VPP(k) = VL(k)+VNL(k). Assuming n electrons per
el
unitcell,theHFdeterminantincludesthe N n orbitalswith
k el III. PERIODICEOMCOUPLED-CLUSTER
thelowesteigenvaluesoutofallk-points. Thetotalenergyis
theusualHFoneplusthenuclear-nuclearrepulsion,whichis
Givenaclosed-shellreferencedeterminant Φ ,theground-
computedusingtheEwaldexpression | (cid:105)
stateCCSDwavefunctionis Ψ =eT Φ whereT =T +T ,
E = 1(cid:88)(cid:48) ZIZJ erfc(ηR R T ) andT1andT2aresingle-and|do0(cid:105)uble-e|xc(cid:105)itationoperato1rs42–425
NN 2 R R T | I − J − |
IJT | I − J − |
+ 214ΩπG(cid:88)(cid:44)0(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88)I ZIeiG·RI(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)2 e−GG22/4η2 (22) T1 =(cid:88)ai k(cid:88)a(cid:48)kitiakkiaEiakkia, (26a)
− √ηπ(cid:88)ZI2− 2Ωπη2 (cid:88)ZI2, T2 = 12(cid:88)abijkak(cid:88)b(cid:48)kikjtiakkiajbkkjbEiakkiaEbjkkbj. (26b)
I I
4
In this work, we use a HF reference determinant, although ofperiodicity,theEOM-CCSDscalingforallchargedexcita-
the equations below and in the Appendix apply to any de- tionsatagivenk-point isn3 n2 N3 forIP-EOM-CCSDand
occ vir k
terminant. In Eqs. (26) and throughout, the indices i, j,k,l n n4 N3 for EA-EOM-CCSD (if the charged excitations at
occ vir k
denote occupied orbitals and a,b,c,d denote virtual orbitals. all sampled k-points are desired, then N such calculations
k
Because the crystal Hamiltonian has translational symmetry, can be performed independently). This specific scaling with
theexcitationoperatorsmustconservecrystalmomentum,i.e. thenumberofoccupiedandvirtualorbitalsusuallymakesthe
(cid:80) (cid:80)
aka− iki = Gwhereka andki arethecrystalmomenta calculation of conduction (virtual) bands significantly more
ofparticleandholeorbitalsandGisareciprocallatticevec- expensivethanthatofvalence(occupied)bands.
tor.Thisrequirementisindicatedbytheprimedsummationin In this work we focus on the quasiparticle excitations.
Eq.(26)andweemphasizethateachprimedsummationindi- These are the many-body states that have the largest over-
catesthatoneofthelistedmomentaisfixedandneednotbe lap with the mean-field single-particle excitations, i.e. the
summed. Introducingthenon-HermitianCCeffectiveHamil- largest r 2 and ra,k2 (a proxy for the EOM Green’s func-
i,k
tonianH¯ e−THeT,theCCSDground-stateenergyandexci- tionpol|est|rength)|. Th|isobservationleadstoanefficienttar-
≡
tationamplitudesaredeterminedby getedDavidsondiagonalizationprocedurebasedonthewave-
functioncharacter,ratherthantheenergy.49,50Ateachk-point
E = Φ H¯ Φ , (27)
0 (cid:104) 0| | 0(cid:105) sampled in the Brillouin zone, we typically target the three
0= ΦakaαH¯ Φ , (28a) lowest-lying IP and EA excitations with such single-particle
(cid:104) ikiα | | 0(cid:105) character. Awayfromthebandedge,thequasiparticlepicture
0=(cid:104)Φaikkiaαα,,jbkkjβbβ|H¯|Φ0(cid:105), (28b) breaks down: the excitations develop an effective linewidth
and nontrivial satellite structure through a growing contribu-
where Φakaα and Φakaα,bkbβ are Slater determinants with one tionofmultipledeterminants,whichalsoleadstoslowercon-
ikiα ikiα,jkjβ
and two electron-hole pairs. The explicit forms of the en- vergenceoftheDavidsonprocedure. Inthisregime,itismore
ergy and amplitude equations with translational symmetry efficienttodirectlyconstructthemomentum-andfrequency-
are given in Ref. 18. Importantly, the computational cost of dependent, one-particle Green’s function as done in our pre-
CCSD for periodic systems scales as n2 n4 N4, where n viousworkontheuniformelectrongas.33
occ vir k occ
and n denote the number of occupied and virtual orbitals
vir
per unit cell; this is a factor of N2 less than the equivalent
k
calculationthatneglectsmomentumconservation. IV. CONVERGENCEANDFINITE-SIZEEFFECTS
Coupled-clusterexcitedstatesandenergiesaredetermined
throughtheequation-of-motion(EOM)formalism,46–48which Itisworthbrieflydiscussingtherequirementsforconverg-
amounts to diagonalizing the effective Hamiltonian H¯ in an ingthecorrelatedelectronicstructureofaperiodicsystem.As
appropriate space of excitations. For a CCSD ground state, inamolecularelectronicstructurecalculation, wemustcon-
wecalculateionizationpotentials(IPs)viadiagonalizationin vergewithrespecttothesingle-particlebasisset,i.e.thesetof
the space of 1-hole (1h) and 2-hole, 1-particle (2h1p) states, Gaussianfunctionsφ˜ (r)inEq.(1),aswellaswithrespectto
µ
andwecalculateelectronaffinities(EAs)viadiagonalization thecorrelationlevel. Convergenceinthesetworegardscanbe
in the space of 1-particle (1p) and 2-particle, 1-hole (2p1h) expectedtobesimilartothatinamolecularsystem.However,
states: |ΨnN,k±1(cid:105) = R±(n,k)|Ψ0(cid:105) = [R±1(n,k) + R±2(n,k)]|Ψ0(cid:105) in addition, we must converge with respect to the Brillouin
whereR (R+)creates1h(1p)excitationsandR (R+)creates zonesampling,i.e.thesetofcrystalmomentumvectorskap-
−1 1 −2 2
2h1p (2p1h) excitations in the IP-EOM (EA-EOM) frame- pearing in Eq. (1). (This may formally be thought of as part
work. Explicitly, ofthebasissetconvergencefortheinfinitesystem). Thereis
substantial experience with this convergence behavior at the
R−α(n,k)=(cid:88)rikaikα+(cid:88) (cid:88)(cid:48) rikibjkkbjEbjkkbjaikiα (29) mreleaatne-dfileelvdelle.2v1e–2l3(,5H1FHaonwdevDeFr,Ti)t,isbuimt mpourctahnltetsosnsooteatththaetccoonr--
i bij kbkikj vergencewithrespecttoBrillouinzonesamplingcanbequite
R+α(n,k)=(cid:88)raka†akα+(cid:88) (cid:88)(cid:48) rakabjkkbja†akaαEbjkkbj (30) sloIwn.general, if all quantities are smooth in the Brillouin
a abj kakbkj zone, then approximating integrals by finite Brillouin zone
sampling has exponentially small error in the grid spacing,
whereagaintheprimedsummationsarerestrictedtoenforce
momentum conservation. The r-amplitudes satisfy the usual i.e. exp( aN−1/3). However, omitting the G = 0 contribu-
− k
IP/EAEOMequations,modifiedtoincludecrystalmomenta; tion in the definition of Coulomb-based integrals is an addi-
explicit expressions are given in the Appendix. Again we tional source of finite-size error. For example, the exact HF
emphasize that each primed summation leads to one crystal exchangeenergyperunitcellisgivenby
momentum which can be fixed by momentum conservation; (cid:90) (cid:90)
(fAor.1ebx)amanpdle(,Ai.n2bth),eofinnlyaltwteormofothfethfoeuRr2craymstpallitmudoemeeqnutaatnioeends EX =−(cid:88)ij BZ Ω(2dπk)3i BZ Ω(2dπk)3j(ikijkj|jkjiki). (31)
to be explicitly summed. Note that although ground-state
CCSDscalesasN6,thesubsequentEOM-CCSDcalculations Inthislimitofinfinitek-pointsampling,thevolumeelement
have a reduced scaling of N5. In particular, in the presence associated with Brillouin zone integration cancels the diver-
5
gentCoulombterm, i.e.itisanintegrabledivergence. How- V. GROUND-STATERESULTS
ever, using the molecular orbital integrals from Eq. (24) and
afiniteBrillouinzonesamplingleadstotheapproximateex-
We first consider the ground-state properties of diamond
changeenergy
and silicon, to establish the convergence properties of the
groundstateCCbeforeproceedingtoourexcitedstatestudies
EX =−NΩk2 (cid:88)ij k(cid:88)ik(cid:48)jG(cid:12)(cid:12)(cid:12)ρiki,jkj(G)(cid:12)(cid:12)(cid:12)2 |G−4πqij|2, (32) itnaltshteruncetxutres,ecwtiiothn.twDoiaamtoomndspanerd(spirliicmointivshe)aruentihteceslalm. Eexccryeps-t
whereotherwisenoted,allcalculationsareperformedwiththe
where the primed summation excludes G = 0 when q = T = 300 K experimental lattice constants a = 3.567 Å and
ij
k k =0. ThisneglectedtermrepresentsanO(1)integrand a=5.431Å,fordiamondandsiliconrespectively.
i j
com−ingfromthenonzerochargedensityρ (G = 0),lead- We use GTH pseudopotentials,40,41 explicitly treating four
ik,ik
ingtoanintegrationerrorofO(N−1/3). Inthesecasesrelated valenceelectronsperatom(eightperunitcell),andmatching
k
toHFexchange,variouscorrectionsexistthataimtoacceler- single-particle basis sets, obtained from the CP2K software
atetheBrillouinzoneconvergence,includingauxiliaryfunc- package.38,39,56 Inthepresentwork, weusepseudopotentials
tion techniques52 and real-space truncation of the Coulomb and basis sets that were originally optimized for use in DFT
interaction.53,54 In this work, all correlated calculations and calculationswiththelocaldensityapproximation(LDA);fu-
their underlying HF calculations are done with no such cor- ture work will consider pseudopotentials and basis sets that
rections. However,whentotalenergiesarerequired,wesep- are optimized for HF calculations24 and we will present a
aratelyconvergetheHFenergyusinganexchangeinteraction comparison with all-electron calculations. The DZV, DZVP,
with a spherically truncated Coulomb potential; these calcu- and TZVP basis sets have 8, 13, and 17 orbitals per atom
lations are performed with isotropic k-point meshes and the (twice as many per cell). The real-space grid spacing used
resultsareextrapolatedusingafinite-sizescalingoftheform fortheintegralswasabout0.17Å,correspondingtoakinetic
exp( aN1/3).53,54 energy cutoff of approximately 400 Ry (note that the same
− k gridisusedtoresolvetheorbitals, Eq.(2), andthepairden-
Finite-sizeerrorsinthecorrelatedtheoriescanbeanalyzed
sities,Eq.(17)). Allcalculationsareconverged,withrespect
inasimilarmanner. Withk-pointsampling, theMP2energy
to this grid, to an accuracy better than 10 4 au per atom. In
perunitcellisapproximatedby −
allground-statecalculations,theBrillouinzonewassampled
E = Ω (cid:88)(cid:48) (cid:88)takabkb ffroormmeaduwniitfhotrhmePΓy-cSeCnFterseodftwgrairde.pAaclklacgael.c5u7lations were per-
MP2 N3 ikijkj
k kakbkikjabij As described above, total convergence in periodic solids
(cid:34) (cid:88)(cid:48) 4π mustbeachievedwithrespecttothesingle-particlebasisand
× 2 ρiki,aka(G) G q 2ρjkj,bkb(Gijab−G) (33) the sampling of the Brillouin zone (finite-size effects). In
G | − ia| Fig. 1 we show the convergence of the ground-state corre-
(cid:35)
(cid:88)(cid:48) 4π lation energy of diamond, at the MP2 and CCSD levels of
− ρiki,bkb(G) G q 2ρjkj,aka(Gijba−G) theory, as a function of the size of the single-particle basis
G | − ib| andthenumberofk-pointssampledintheBrillouinzone. As
expected from molecular calculations with comparable basis
whereagainthetermwithG =0isneglectedwhenq =k
ia i− sets, the totalcorrelationenergies arenotconverged withre-
k =0(firstterm)orq =k k =0(secondterm). Unlike
a ib i− b specttothesingle-particlebasis;thecorrelationenergydiffer-
in the case of exchange, these neglected terms vanish at the
encebetweenresultsobtainedwiththeDZVPandTZVPbasis
originq = 0becausetheorthogonalityoforbitalsguarantees
that ρ (G = 0) = 0; instead, these terms are associated
ik,ak
with an integrand of O(q2), leading to an integrated error of MP2 CCSD
O(N 1). If the t-amplitudes are correct, then this is the only --44 --44
k−
integrationerror. However,theamplitudes ell) --55 --55
c
/V
takabkb = (ikiaka|jkjbkb)∗ (34) y(e --66 --66
ikijkj εiki +εjkj −εaka −εbkb nerg --77 --77
e
n
inherittheerroroftheHForbitalenergies,whichcanbeana- atio --88 --88
lyzedanalogously. Orthogonalityguaranteesthatunoccupied el DZV
energieshavethesamefavorable Nk−1 error; however, justas Corr --99 --99 DTZZVVPP
fortheexchangeenergy,theoccupiedorbitalsexhibitanN−1/3 --1100 --1100
k 11 22 33 44 11 22 33 44
error.55 Without any corrections to the orbital energies, this (Numberofk-points)1/3 (Numberofk-points)1/3
pollutes the t-amplitudes and thus dominates the error in the
correlation energies (see below). The inclusion of all four- FIG.1.CorrelationenergyperunitcellofdiamondcalculatedwithMP2(left,
ienrrdoerx. integrals in CCSD is similarly responsible for a Nk−1/3 oNpke→ns∞ymebxotrlsa)poalnadteCdCvaSlDue(sraigshsut,mfiilnlegdasfiynmitbeo-slisz)e. sDcaalsihnegdolfintheesfionrdmicaNtek−1t/h3e.
6
DZV DZVP TZVP a(Å) B(GPa) ∆E(eV/atom)
HF C
Extrap,exp( aN1/3) -301.05 -301.86 -301.95 HF 3.527 (3.552) 507 (495) 5.36 (5.28)
− k
MP2 3.545 (3.553) 436 (450) 7.91 (7.97,8.039)
MP2 CCSD 3.539 — 463 — 7.04 (7.295)
3 3 3 -5.91 -8.27 -8.81 Experiment 3.553 455 7.55
× ×
4 4 4 -5.96 -8.24 -8.70
Ex×trap×,Nk−1/3 -6.10 -8.16 -8.36 Si
HF 5.435 (5.512) 107 (103) 3.03 (2.97)
CCSD MP2 5.347 (5.415) 101 (100) 4.96 (5.05)
3 3 3 -5.01 -7.00 -7.45 CCSD 5.393 — 103 — 4.15 —
× ×
4 4 4 -5.15 -7.12 -7.55 Experiment 5.421 101 4.68
Ex×trap×,Nk−1/3 -5.58 -7.50 -7.86
TABLEII.Latticeconstanta,bulkmodulusB,andcohesiveenergy∆E of
TABLEI.Hartree-FocktotalenergyperunitcellandMP2/CCSDcorrelation diamondandsiliconusingtheTZVPsingle-particlebasis. Latticeconstant
energyperunitcellofdiamondineV. andbulkmoduluswereobtainedfromfitstocalculationswitha3×3×3k-
pointmesh.Crystalenergycontributionstothecohesiveenergyarecalculated
attheT =300KlatticeconstantandextrapolatedtoNk →∞.Valuesgiven
in parentheses are those calculated with the plane-wave PAW approach as
setsis0.4–0.5eVperunitcell(halfthatperatom),tobecom- reportedinRefs.22and25.Allexperimentalvalueshavebeencorrectedfor
paredtotheHFenergydifferences,whicharelessthan0.1eV zero-pointvibrationaleffects.59
(notshown). Inthiswork,wedonotpursueacompletebasis
set limit extrapolation based only on double-zeta and triple-
zetabasissets. Itcanbeexpected,however,thatthisbasisset 7.295eV/atom.25 Similarly,ourMP2cohesiveenergyforsil-
errorislocalincharacterandcanbecorrectedusingstandard icon, 4.96 eV/atom, is in good agreement with the previ-
explicitr correlationtechniques.58Furthermore,weobserve ouslyreportedMP2value,5.05eV/atom.22Thediscrepancies
12
thatthebasissetcorrectionstothecorrelationenergyaresim- in our data of order 0.1–0.2 eV/atom likely originate from
ilarforMP2andCCSD(asinmolecules),andwecouldthus a combination of the finite basis set and treatment of core
similarlyusecompletebasissetlimitMP2calculationstocor- electrons. The calculations reported in Refs. 22 and 25 em-
recttheCCSDresults. ployaplane-wavebasissetandtheprojector-augmentedwave
Withregardstofinite-sizeeffects,weobservethatinsmall (PAW)method,whichdiffersfromthenorm-conservingpseu-
dopotentialsusedhere. Furthermore,wealwaysfindthatthe
basissetstheMP2correlationenergyconvergesmorequickly
absolutevalueofthecorrelationenergygivenbyMP2islarger
than the CCSD energy with the number of k-points. How-
than that given by CCSD, which is the opposite of that ob-
ever,forthelargest(TZVP)basisset,theenergyconvergence
served in Ref. 25; this difference may also be attributable to
issimilar. InTab.I,wegivetheHFandcorrelationenergies
thedifferencebetweenthePAWandpseudopotentialapprox-
obtainedwitha3 3 3and4 4 4samplingoftheBril-
× × × × imations.
louinzone;thelargestCCcalculationpresentedhere,usinga
4 4 4samplingoftheBrillouinzoneandtheTZVPbasis, Usinga3 3 3k-pointmeshandtheTZVPbasis,wehave
co×nsti×tutesacanonicalCCSDcalculationof256electronsin furthercalcu×late×dtheequationofstateofdiamondandsilicon,
2,176 orbitals, demonstrating the savings provided by incor- showninFig.3.Weperformedatleastsevenground-statecal-
poratingperiodicsymmetry.
InTab.I,wealsoreportthethermodynamiclimitbasedon
Diamond Silicon
an extrapolation of the form N−1/3, as discussed in Sec. IV. 88..44 55..22
This extrapolation is graphicaklly demonstrated in Fig. 2, m) 88..22 55..00
whichshowsthecohesiveenergy(theenergy-per-atomdiffer- ato 88..00 44..88
ence between the crystal and isolated atom) of diamond and /V 77..88 44..66
(e 77..66
siliconasafunctionofthenumberofk-pointssampledusing y 44..44
the TZVP basis. Specifically, the fitting is restricted to the energ 7777....2424 44..22
tAuwntooremlsatirrcigceetnsetderisgMoietPrso2opaficnodkp-eCpnoC-isSnhDtemllienCshtahenesd,finSeil×dwnoefr×ecncryawslctiatuhlllainntee=dbw3a,ist4ihs. Cohesive 667667......680680 334334......680680 MCCPS2D
functions (to account for basis set superposition error) and 66..44 33..44
the HF contribution to the crystal energy is separately con- 44−−11 33−−11 22−−11 44−−11 33−−11 22−−11
vergedwithrespecttothenumberofk-points, inthemanner (Numberofk-points)−1/3 (Numberofk-points)−1/3
describedinSec.IV.Predictedcohesiveenergiesaregivenin
FIG.2.Cohesiveenergyofdiamond(left)andsilicon(right),calculatedwith
Tab.II.
theTZVPbasis.Foreachk-pointmesh,MP2andCCSDcorrelationenergies
Our cohesive energies for diamond, 7.91 eV/atom (MP2) areaddedtoaseparatelyconvergedHFcrystalenergy.Thecorrelationcom-
and7.04eV/atom(CCSD),areinreasonableagreementwith ponentofthecohesiveenergyisextrapolatedassumingafinite-sizeerrorthat
thoseofpreviousperiodicMP2calculations, 7.97eV/atom22 scaleslikeNk−1/3.Zero-pointcorrectedexperimentalresults59areshownasa
and 8.039 eV/atom,25 and of a previous CCSD calculation, thingreyline.PredictedandexperimentalvaluesaregiveninTab.II.
7
Diamond Silicon wecomputethequasiparticleexcitationsintheionizationpo-
00..2244 00..1166
V) MP2 tential(IP)andelectronaffinity(EA)chargedsectors. Toob-
e 00..2211 00..1144
ell( 00..1188 00..1122 CCSD taincontinuousbandsinconventionalbandstructurecalcula-
c tionswithHForDFT,aninitialself-consistentcalculationis
nit 00..1155 00..1100
u performed with a fixed (and potentially coarse) k-point grid,
per 00..1122 00..0088 which provides an approximate representation of the density
gy 00..0099 00..0066 matrix or density; see Eq. (13). From this latter object, the
er
en 00..0066 00..0044 one-bodyFock-likematrixcanbeconstructedatarbitraryk-
e
v 00..0033 00..0022 pointsanddiagonalized,leadingtonon-self-consistenteigen-
ati
el 00..0000 00..0000 values. To construct a continuous band-structure in EOM-
R
CCSD, we here perform a large series of (independent) cal-
1100..001100..551111..001111..551122..00 3366..00 3388..00 4400..00 4422..00
culationswithak-pointmeshthatisshiftedtoincludethede-
Unitcellvolume(Å3) Unitcellvolume(Å3)
sired k-point. Although this means that certain calculations
FIG.3.Equationofstateofdiamond(left)andsilicon(right)withtheTZVP along the band path will be performed with a sub-optimal
basisanda3 3 3samplingoftheBrillouinzone. mesh,theresultswillconvergeproperlyinthelimitofdense
× ×
k-point sampling. We only present results here for a band
structure path through the Brillouin zone, but the same tech-
culationswithunitcellvolumesvaryingbyabout 10%from niquecouldbeappliedtocalculatethefullquasiparticleden-
±
theequilibriumvalueandfittheresultstoathird-orderBirch- sityofstates.
Murnaghanform,leadingtopredicted(zero-temperature)lat- The band structures of diamond and silicon are shown in
tice constants and bulk moduli, given in Tab. II. Again, the Figs. 4 and 5 using DFT with the PBE functional,60 HF, and
HF and MP2 structural predictions are in good agreement EOM-CCSD using the DZVP basis and a 3 3 3 k-point
with previous values,22 and the observed discrepancy at the mesh. Qualitatively, the EOM-CCSD result×repr×oduces the
HF level provides some measure of the error incurred due expectedbehavior:itpredictsabandgapthatisslightlylarger
to the use of pseudopotentials. Ultimately, the overall good than that of PBE and significantly smaller than that of HF.
agreement between our results and those of Refs. 22 and Similarly, the EOM-CCSD band widths are intermediate be-
25 suggests that the correlation energy differences required tweenthetwomean-fieldresults.
forstructuralpropertiesandcohesiveenergiesarereasonably Tostudytheconvergenceoftheexcited-statestructure,we
converged in both implementations, with respect to both the present in Fig. 6 the calculated value of the indirect (mini-
basisandk-pointsampling. Atleastforthetwosemiconduc- mum) band gap in diamond and silicon, as a function of the
tors studied here, MP2 and CCSD provide a similar level of basis set and number of k-points sampled. Since we do not
accuracy when compared to experimental values determined calculate the full band path due to its high cost, the indirect
byground-stateenergetics. bandgapisevaluatedfromtwoindependentcalculations: an
IP-EOM-CCSD calculation with a Γ-centered mesh and an
EA-EOM-CCSDcalculationwithak-pointmeshcenteredat
VI. EXCITED-STATERESULTS 85%oftheΓ Xline(nearthezoneboundary). Weseethat
∼ −
larger single-particle basis sets favor a smaller bandgap and
finer k-point meshes favor a larger band gap. The band gap
We now turn to the excited-state and spectral calculations
difference between calculations using the DZVP and TZVP
using the EOM-CCSD formalism. As discussed in Sec. III,
basis sets is about 0.1 eV, and relatively independent of the
numberofk-pointssampled.
PBE HF EOM-CCSD
2200
PBE HF EOM-CCSD
1155 1122
1100
1100
88
V)
(e 55 66
ergy 00 eV) 44
n (
E gy 22
55 er
−− En 00
−−1100 −−22
−−1155LL ΓΓ XX WW KK ΓΓ −−44
66
−−
FIG. 4. Band structure of diamond calculated with DFT (PBE), HF, and LL ΓΓ XX WW KK ΓΓ
EOM-CCSD,usingtheDZVPsingle-particlebasisanda3 3 3k-point
× ×
mesh. FIG.5.ThesameasinFig.4,butforsilicon.
8
55..66 Diamond 11..33 Silicon VII. DISCUSSIONANDCONCLUSION
55..44 11..22
V) In this work, we have presented the first-ever results from
(e 55..22 11..11 canonical Gaussian-based coupled-cluster theory for three-
p
a
g 55..00 11..00 dimensional solids. We discussed concomitant finite-size er-
d
ban 44..88 00..99 rors in the mean-field and correlated calculations. Based
ct on ground-state coupled-cluster calculations, we presented
dire 44..66 00..88 the lattice constant, bulk modulus, and cohesive energies of
In 44..44 00..77 DZVP diamond and silicon; these results are in good agreement
TZVP
with limited existing data based on the plane-wave PAW ap-
44..22 00..66
44−−11 33−−11 22−−11 44−−11 33−−11 22−−11 proach.22,25 We further presented the first results of excited-
(Numberofk-points)−1/3 (Numberofk-points)−1/3 state(equation-of-motion)coupled-clustertheoryfortheband
structureandbandgapofthesametwosemiconductors. Our
FIG.6.Indirectbandgapofdiamond(left)andsilicon(right)calculatedwith predicted minimum band gaps, 5.37 eV for diamond and
EOM-CCSD.DashedlinesindicateextrapolationoftheformNk−1/3basedon 1.19 eV for silicon, are in good agreement with experimen-
theisotropic3 3 3and4 4 4samplingmeshes;greyhorizontallines
× × × × talvalues.
indicateexperimentalvalues. FortheTZVPdata,theextrapolationpredicts
fordiamondabandgapof5.37eVcomparedtotheexperimentalvalueof Agreatdealofworkremainstobedonetomakecoupled-
5.45–5.50eV,61,62 andforsiliconabandgapof1.19eV,comparedtothe clustercalculationsonweakly-correlatedsolidsroutine.Infu-
experimentalvalueof1.17eV.62,63 turework,wewillcarefullycomparetheresultsofpseudopo-
tential and all-electron calculations, which represents an un-
quantifiedsourceoferror. Paralleleffortsareaimedatreduc-
The band gap convergence with the number of k-points ingthecostandincreasingourunderstandingofthesecalcu-
is slow and oscillatory (when including anisotropic k-point lations,forexamplethroughlocalcorrelationapproaches,ex-
meshes), a behavior inherited from the underlying HF band plicitly correlated formulations, perturbative corrections due
structure (not shown). In particular, we retain the slow con- totripleexcitations,andapplicationstoawiderrangeofinsu-
vergence arising from the G = 0, q 0 contributions to latingandmetallicmaterials.
→
the energy discussed in Sec. II, which appears to lead to an
N−1/3 scaling of the error in the bandgap. Similarly to the
k
ground-stateenergyhowever, thefinite-sizeeffectsaremuch ACKNOWLEDGMENTS
smaller in the correlated band gap than in the HF band gap:
in diamond, the HF band gap changes by about 7 eV during J.M., Q.S., and G.K.-L.C. were supported by the US De-
convergence,whiletheEOM-CCSDbandgapchangesbyless partment of Energy, Office of Science under Award Num-
than1eV. ber DE-SC0008624. Additional support was provided by
Likeinourtreatmentoftheground-stateandcohesiveen- the Simons Foundation, via the Collaboration on the Many-
ergy,wehaveappliedaN−1/3 finite-sizescalingtotheEOM- ElectronProblem,andaSimonsInvestigatorshipinTheoreti-
k
CCSDbandgaps,restrictedtothetwolargestisotropick-point calPhysics. T.C.B.wassupportedbythePrincetonCenterfor
meshes. Within the approximations of our approach (GTH TheoreticalScienceandbystart-upfundsfromtheUniversity
pseudopotentials, finite basis, and EOM-CCSD correlation) ofChicago.
we predict band gaps of 5.37 eV and 1.19 eV, for diamond
and silicon respectively; the experimental T = 0 band gaps
are5.45–5.50eV61,62and1.17eV,62,63whicharecorrectedfor
excitonic and finite-temperature effects, but not for the more
challenging electron-phonon (vibrational) effects. For dia-
mondinparticular, azero-pointbandgaprenormalizationof
0.4–0.6 eV has been estimated,64–66 which would suggest an
experimentalelectronicbandgapof5.8eVormore.Forcom-
parison,tightlyconvergedall-electronG W calculationspre-
0 0
dictbandgapsof5.61–5.62eVand1.10–1.11eV,67,68 which
are quite close to G W calculations with norm-conserving
0 0
pseudopotentials,69 but 0.1–0.2 eV smaller than all-electron
self-consistentGWcalculations70and0.2–0.3eVsmallerthan
vertex-corrected calculations.12 Interestingly, although MP2
predictsground-statestructuralpropertiesthatareonparwith
CCSD,itmassivelyunderestimatesthebandgapsofthesetwo
materials, giving 1.9 eV for diamond and 1.2 eV for sili-
−
con,22whichwasattributedtothelargepolarizabilityofthese
materials.
9
Appendix:IP/EA-EOM-CCSDequationsforperiodicsystems
Assumingaclosed-shellreference,thespatial-orbitalIP-EOM-CCSDamplitudeequationsaregivenby71
(cid:88)(cid:88) (cid:88)(cid:88)
(H¯R−)iki = (cid:48) −Ukkkikirkkk + (cid:48) Ulkldkd(2rikildkkld −rlklidkkid)
k kk ld klkd (A.1a)
(cid:88) (cid:88)
(cid:48) (2W W )r dkd
− kkklklikidkd − lklkkkikidkd kkklkl
kld kkklkd
(cid:88)(cid:88) (cid:88)(cid:88) (cid:88)(cid:88)
(H¯R−)ikibjkkbj = (cid:48) Ubkbdkdrikidjkkjd − (cid:48) Ulklikirlklbjkkbj − (cid:48) Ulkljkjrikilbkklb
d kd l kl l kl
(cid:88)(cid:88) (cid:88)(cid:88) (cid:88)(cid:88)
(cid:48) W r + (cid:48) (2W W )r dkd (cid:48) W r dkd
− kkkbkbikijkj kkk lklbkbdkdjkj − bkblkldkdjkj ikilkl − lklbkbdkdjkj lkliki
k kk ld klkd ld klkd (A.1b)
(cid:88)(cid:88) (cid:88)(cid:88)
+ (cid:48) W r bkb (cid:48) W r dkd
kkklklikijkj kkklkl − kkkbkbikidkd kkkjkj
kl kkkl kd kkkd
 
−(cid:88)(cid:88)(cid:48) (cid:88) (cid:88)(cid:48) (2Wlklkkkdkdckc −Wkkklkldkdckc)rkkkldkkldtickkicjbkkjb,
c kc kld kkklkd
andtheEA-EOM-CCSDequationsaregivenby72
(cid:88)(cid:88) (cid:88)(cid:88)
(H¯R+)aka = (cid:48) Uakackcrckc + (cid:48) Ulkldkd(2rakaldkkld −rdkdlakkla)
c kc ld klkd (A.2a)
(cid:88) (cid:88)
+ (cid:48) (2W W )rckcdkd
akalklckcdkd − akalkldkdckc lkl
lcd klkckd
(cid:88)(cid:88) (cid:88)(cid:88) (cid:88)(cid:88)
(H¯R+)akabkb = (cid:48) U rakabkb + (cid:48) U rckcbkb + (cid:48) U rakadkd
jkj − lkljkj lkl akackc jkj bkbdkd jkj
l kl c kc d kd
(cid:88)(cid:88) (cid:88)(cid:88) (cid:88)(cid:88)
+ (cid:48) Wakabkbckcjkjrckc + (cid:48) (2Wlklbkbdkdjkj −Wbkblkldkdjkj)rakaldkkld − (cid:48) Wakalklckcjkjrckclbkklb
c kc ld klkd lc klkc
(cid:88)(cid:88) (cid:88)(cid:88)
(cid:48) W rckcaka + (cid:48) W rckcdkd
− bkblkljkjckc lkl akabkbckcdkd jkj
lc klkc cd kckd
 
−(cid:88)(cid:88)(cid:48) (cid:88) (cid:88)(cid:48) (2Wkkklklckcdkd −Wkkklkldkdckc)rckcldkkldtkakkkabjkkjb.
k kk lcd klkckd
(A.2b)
Asdescribedinthetext,primedsummationsindicatethatoneofthelistedmomentaisfixedbymomentumconservation. The
intermediatesUandWaretheusualonesarisingincoupled-clustertheoryfromtheeffectiveHamiltonianH¯,73andtheirefficient
calculationalsoaccountsformomentumconservation.
1J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 8L. Reining, V. Olevano, A. Rubio, and G. Onida, Phys. Rev. Lett. 88,
(2003). 066404(2002).
2S.J.ClarkandJ.Robertson,Phys.Rev.B82,085208(2010). 9M.RohlfingandS.G.Louie,Phys.Rev.B62,4927(2000).
3J.M.Crowley,J.Tahir-Kheli, andW.A.Goddard,J.Phys.Chem.Lett.7, 10M.Shishkin, M.Marsman, andG.Kresse,Phys.Rev.Lett.99,246403
1198(2016). (2007).
4F.TranandP.Blaha,Phys.Rev.Lett.102,226401(2009). 11P. Romaniello, S. Guyot, and L. Reining, J. Chem. Phys. 131, 154111
5M.S.HybertsenandS.G.Louie,Phys.Rev.Lett.55,1418(1985). (2009).
6M.S.HybertsenandS.G.Louie,Phys.Rev.B34,5390(1986). 12A. Gru¨neis, G. Kresse, Y. Hinuma, and F. Oba, Phys. Rev. Lett. 112,
7P.L.deBoeij,F.Kootstra,J.A.Berger,R.vanLeeuwen, andJ.G.Snijders, 096401(2014).
J.Chem.Phys.115,1995(2001). 13Y.Pavlyukh,A.-M.Uimonen,G.Stefanucci, andR.vanLeeuwen,Phys.
Rev.Lett.117,206402(2016).
10
14S.Hirata,M.Head-Gordon, andR.J.Bartlett,J.Chem.Phys.111,10774 41C.Hartwigsen,S.Goedecker, andJ.Hutter,Phys.Rev.B58,3641(1998).
(1999). 42J.Cˇ´ızˇek,J.Chem.Phys.45,4256(1966).
15J.SunandR.J.Bartlett,J.Chem.Phys.104,8553(1996). 43G.E.Scuseria,C.L.Janssen, andH.F.SchaeferIII,J.Chem.Phys.89,
16P.Y.Ayala,K.N.Kudin, andG.E.Scuseria,J.Chem.Phys.115,9698 7382(1988).
(2001). 44I.ShavittandR.J.Bartlett,Many-BodyMethodsinChemistryandPhysics:
17S.Hirata,I.Grabowski,M.Tobita, andR.J.Bartlett,Chem.Phys.Lett. MBPT and Coupled-Cluster Theory (Cambridge University Press, Cam-
345,475 (2001). bridge;NewYork,2009).
18S.Hirata,R.Podeszwa,M.Tobita, andR.J.Bartlett,J.Chem.Phys.120, 45R.J.BartlettandM.Musiał,Rev.Mod.Phys.79,291(2007).
2581(2004). 46H.J.Monkhorst,Int.J.QuantumChem.12,421(1977).
19C.Pisani,J.Mol.Struct.(Theochem)621,141(2003). 47J.F.StantonandR.J.Bartlett,J.Chem.Phys.98,7029(1993).
20C. Pisani, M. Busso, G. Capecchi, S. Casassa, R. Dovesi, L. Maschio, 48A.I.Krylov,Ann.Rev.Phys.Chem.59,433(2008).
C.Zicovich-Wilson, andM.Schu¨tz,J.Chem.Phys.122,094113(2005). 49W.ButscherandW.Kammer,J.Comput.Phys.20,313(1976).
21M.Marsman, A.Gru¨neis, J.Paier, andG.Kresse,J.Chem.Phys.130, 50D.Zuev,E.Vecharynski,C.Yang,N.Orms, andA.I.Krylov,J.Comp.
184103(2009). Chem.36,273(2015).
22A. Gru¨neis, M. Marsman, and G. Kresse, J. Chem. Phys. 133, 074107 51K.LiaoandA.Gru¨neis,J.Chem.Phys.145,141102(2016).
(2010). 52F.GygiandA.Baldereschi,Phys.Rev.B34,4405(1986).
23A.Gru¨neis,G.H.Booth,M.Marsman,J.Spencer,A.Alavi, andG.Kresse, 53J.SpencerandA.Alavi,Phys.Rev.B77,193110(2008).
J.Chem.TheoryComput.7,2780(2011). 54R.SundararamanandT.A.Arias,Phys.Rev.B87,165122(2013).
24M.DelBen,J.Hutter, andJ.VandeVondele,J.Chem.TheoryComput.8, 55G.AissingandH.J.Monkhorst,Int.J.Quant.Chem.48,81(1993).
4177(2012). 56J.VandeVondeleandJ.Hutter,J.Chem.Phys.127,114105(2007).
25G.Booth,A.Gru¨neis,G.Kresse, andA.Alavi,Nature493,365(2013). 57“PySCF version 1.2,” (2016), PySCF is freely available from
26H.Stoll,Phys.Rev.B46,6700(1992). https://github.com/sunqm/pyscf.
27H.Stoll,B.Paulus, andP.Fulde,J.Chem.Phys.123,144108(2005). 58W.Klopper, F.R.Manby, S.Ten-No, andE.F.Valeev,Int.Rev.Phys.
28B.Paulus,Phys.Rep.428,1 (2006). Chem.25,427(2006).
29S.J.Nolan,M.J.Gillan,D.Alfe`,N.L.Allan, andF.R.Manby,Phys.Rev. 59L.Schimka,J.Harl, andG.Kresse,J.Chem.Phys.134,024116(2011).
B80,165109(2009). 60J.P.Perdew,K.Burke, andM.Ernzerhof,Phys.Rev.Lett.77,3865(1996).
30F.Libisch,C.Huang, andE.A.Carter,Acc.Chem.Res.47,2768(2014). 61C.D.Clark,P.J.Dean, andP.V.Harris,Proc.R.Soc.Lond.A277,312
31G.Kotliar,S.Y.Savrasov,K.Haule,V.S.Oudovenko,O.Parcollet, and (1964).
C.A.Marianetti,Rev.Mod.Phys.78,865(2006). 62K.P.O’DonnellandX.Chen,Appl.Phys.Lett.58,2924(1991).
32Q.SunandG.K.-L.Chan,Acc.Chem.Res.49,2705(2016). 63W.Bludau,A.Onton, andW.Heinke,J.Appl.Phys.45,1846(1974).
33J.McClain,J.Lischner,T.Watson,D.A.Matthews,E.Ronca,S.G.Louie, 64M.Cardona,SolidStateComm.133,3(2005).
T.C.Berkelbach, andG.K.-L.Chan,Phys.Rev.B93,235139(2016). 65F.Giustino,S.G.Louie, andM.L.Cohen,Phys.Rev.Lett.105,265501
34L.Hedin,Phys.Scr.21,477(1980). (2010).
35F. Aryasetiawan, L. Hedin, and K. Karlsson, Phys. Rev. Lett. 77, 2268 66G.Antonius,S.Ponce´,P.Boulanger,M.Coˆte´, andX.Gonze,Phys.Rev.
(1996). Lett.112,215501(2014).
36L.Hedin, B.I.Lundqvist, andS.Lundqvist,SolidStateComm.5,237 67C. Friedrich, S. Blu¨gel, and A. Schindlmayr, Phys. Rev. B 81, 125102
(1967). (2010).
37B.I.Lundqvist,Phys.Kondens.Mater.6,193(1967). 68J.Klimesˇ,M.Kaltak, andG.Kresse,Phys.Rev.B90,075125(2014).
38J.VandeVondele,M.Krack,F.Mohamed,M.Parrinello,T.Chassaing, and 69W.ChenandA.Pasquarello,Phys.Rev.B92,041115(2015).
J.Hutter,Comput.Phys.Commun.167,103 (2005). 70M.ShishkinandG.Kresse,Phys.Rev.B75,235102(2007).
39TheCP2Kdevelopersgroup,“CP2Kversion3.0,” (2016),CP2Kisfreely 71M.NooijenandJ.G.Snijders,Int.J.QuantumChem.48,15(1993).
availablefromhttp://cp2k.org/. 72M.NooijenandR.J.Bartlett,J.Chem.Phys.102,3629(1995).
40S.Goedecker,M.Teter, andJ.Hutter,Phys.Rev.B54,1703(1996). 73M.NooijenandJ.G.Snijders,Int.J.QuantumChem.47,3(1993).