The elphbolt ab initio solver for the coupled electron-phonon Boltzmann transport equations

ARTICLE OPEN The elphbolt ab initio solver for the coupled electron-phonon Boltzmann transport equations Nakib H. Protik 1 ✉, Chunhua Li 2 , Miguel Pruneda 1 , David Broido 2 and Pablo Ordejón 1 ✉ elphbolt is a modern Fortran (2018 standard) code for efficiently solving the coupled electron–phonon Boltzmann transport equations from first principles. Using results from density functional and density functional perturbation theory as inputs, it can calculate the effect of the non-equilibrium phonons on the electronic transport (phonon drag) and non-equilibrium electrons on the phononic transport (electron drag) in a fully self-consistent manner and obeying the constraints mandated by thermodynamics. It can calculate the lattice, charge, and thermoelectric transport coefficients for the temperature gradient and electric fields, and the effect of the mutual electron–phonon drag on these transport properties. The code fully exploits the symmetries of the crystal and the transport-active window to allow the sampling of extremely fine electron and phonon wave vector meshes required for accurately capturing the drag phenomena. The coarray feature of modern Fortran, which offers native and convenient support for parallelization, is utilized. The code is compact, readable, well-documented, and extensible by design. npj Computational Materials (2022)8:28 ; https://doi.org/10.1038/s41524-022-00710-0 INTRODUCTION Ab initio computation of the transport properties of materials allows a deeper understanding and probing of the rich underlying physics. It is also crucial for the predictive designing of materials for industrial applications. The advances in density functional theory (DFT) 1,2 and density functional perturbation theory (DFPT) 3 have enabled accurate electron (e) and phonon (ph) band structure calculations. Furthermore, various freely available codes exist that allow fast computation of the ph–ph 4 and e–ph 5,6 interactions. On the transport side, significant progress has been made that allows computations of the mode (band/branch and wave vector) resolved lifetimes of electrons and phonons. The phonon thermal and the electronic charge conductivity can now be calculated ab initio in the relaxation time approximation (RTA) or, going one step further, from a full solution of the corresponding single-species Boltzmann transport equation (BTE) 4,6–9 . In an interacting e–ph gas, however, the transport of the two systems is intimately connected, and a unified, coupled e–ph transport may, instead, be needed for an accurate description of the underlying physics. The idea of a complete separation of the e and ph BTEs dates back to the 1930s and is known as Bloch’s Assumption 10 .A paradoxical consequence of this assumption is that when one solves the e (ph) BTE, the phonon (electron) system is taken to remain in equilibrium. This framework was promptly questioned by Peierls 11 , who argued that there must exist a momentummixing between the electrons and the phonons, which would, in turn, cause both systems to move under the influence of a driving field. A theory of the coupled e–ph transport was developed by Gurevich in 1946 12 , which allowed calculation of the effect of the non-equilibrium electrons (phonons) on the phonons (electrons). The first is called the electron drag and the latter, phonon drag. In an interacting electron–phonon gas, however, the two effects are inseparable and a mutual electron–phonon drag effect occurs. A decade later the first experimental evidence for the drag effect on the thermopower of germanium and silicon was found by Frederikse 13 and Geballe and Hull 14,15 . Around the same time, an influential theory was developed by Conyers Herring 16 to explain the experimental findings. Herring’s theory, which contains several free-parameters, partially decouples the e and the ph BTEs, retaining an approximate term to account for the drag effect. From a physical point of view, however, Herring’s theory violates a deep, thermodynamic connection between the Seebeck and the Peltier effects known as the Kelvin–Onsager relationship 17 . Semi-analytical work on 2D systems was carried out in 1987 by Cantrell, Butcher, and coworkers 18,19 . More recently, in 2014, a semi-analytical model with ab initio fitted parameters and a partially decoupled framework was employed by Mahan, Lindsay, and Broido to calculate the drag effect on the thermopower of silicon 20 . Soon after, fully ab initio drag calculations with partially decoupled solutions of the e and ph BTEs were devised by Zhou et. al. in 2015 for silicon 21 . The code used from that work was released in 2020 22 . Ab initio solutions to partially decoupled e and ph BTEs were also developed by Fiorentini and Bonini in 2016 for silicon 23 and Macheda and Bonini 24 in 2018 for diamond. Finally, in 2020, a fully coupled e–ph BTEs solution was devised and applied to gallium arsenide using model e–ph interactions by Protik and Broido 25 . In the same year, that method was improved to include fully ab initio e–ph interactions to calculate the drag effect in silicon carbide by Protik and Kozinsky 26 . Here we present elphbolt (short for electron–phonon Boltzmann transport), a code that features major improvements over the methods given in refs. 25,26 . Moreover, we release the code as Free/Libre software under the GNU General Public License version 3, bringing the capabilities for calculating the e–ph drag physics via an ab initio solution of the fully coupled e–ph BTEs, almost a century after Peierls’conception of the idea, to the broader transport physics community. Our code is hosted on github 27 . Using an ab initio and Kelvin–Onsager relationship conserving solution of the coupled e–ph BTEs, elphbolt gives access to the: 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus Bellaterra, 8193 Barcelona, Spain. 2 Department of Physics, Boston College, Chestnut Hill, Boston, MA 02467, USA. ✉email: [email protected]; pablo.or[email protected] www.nature.com/npjcompumats Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences 1234567890():,; ●mode resolved phonon thermal conductivity; ●mode resolved phonon Peltier coefficient; ●mode resolved electronic charge conductivity; ●mode resolved electronic thermal conductivity; ●mode resolved electronic Seebeck and Peltier coefficients; ●and the effect of e–ph drag on all of the above quantities. This code is suitable for the study of 3dand 2dinsulators, semiconductors, semimetals, and metals. In the sections below, we present the theory behind elphbolt, the implementation, and the outlook. RESULTS In this section, we give results for the calculated thermopower, mobility, and thermal conductivity of n-doped silicon. First, the basic ingredients of the theory are described in detail. We present the elementary interactions considered in this work. This is followed by a discussion of the BTEs in the forms in which they are implemented, and the various types of solutions that elphbolt offers. Lastly, we present the transport coefficients that are obtained from the solution of the BTEs along with a brief discussion of the Kelvin–Onsager reciprocal relationship connecting the thermoelectric coefficients. Electrons, phonons, and interactions First, we need the electrons and phonons, and their interactions calculated on arbitrarily fine wave vector meshes. These are achieved by the Wannier interpolation techniques described in refs. 5,28–30 . We refer the readers to these seminal works for the details of the calculation of the Wannier functions and obtaining the Wannier representations of the various physical quantities starting from their Bloch representations. In particular, the expressions for the Hamiltonian, dynamical matrix, and the e–ph matrix elements are given in refs. 5,30 . Below we simply provide the expressions for the Wannier to Bloch transformations that are computed within elphbolt. The Hamiltonian in the Bloch representation can be obtained by the Fourier transformation: HmnðkÞ¼ 1 NeX Re expðikReÞHmnðReÞ;(1) where mand nare band indices, kis an arbitrary wave vector, H(R e ) is the Hamiltonian in the Wannier representation, living on the real space spanned by {R e }, and N e is the number of real space unit cells. Diagonalizing H(k), we obtain the electronic band energies ϵ mk . The unitary matrix U k digonalizing the Hamiltonian contains the eigenstates mk ji . The band velocities are obtained from the Hellmann–Feynman theorem 31 : vmk¼1 _mkhj∇kHðkÞmkji;(2) where ℏis the reduced Planck constant. Similarly, the dynamical matrix in the Bloch representation is obtained from its Wannier representation using Dss0ðqÞ¼ 1 Nph X Rph expðiqRphÞDss0ðRphÞþDNAC ss0ðqÞ;(3) where sand s0denote the phonon branches, qis an arbitrary wave vector, Dss0ðRphÞis the dynamical matrix in the Wannier representation, and R ph locates each unit cell in real space containing N ph cells. The first term is short-ranged, and the term DNAC ss0is the long-range (the, so-called, nonanalytic) correction due to the dipole–dipole interaction given by 32 DNAC ss0ðqÞ¼ 1 ffiffiffiffiffiffiffiffiffiffiffiffiffi mτmτ0 pe2 ε0V qZ τqZ τ0 qϵ1q;(4) where eis the electronic charge, Vis the primitive unit cell volume, τand τ0label the basis atoms with masses m τ and mτ0,Z * is the Born effective charge tensor, ε 0 is the permittivity of free space, and ϵ ∞ is the high-frequency dielectric tensor. This additional term is only required for polar materials. We obtain the phonon branch energies ℏω sq and the eigenstates sq ji by diagonalizing Dss0ðqÞwith the unitary matrix u q . Here ω sq is the phonon angular frequency. The phonon group velocities are calculated using vsq¼1 2ωsq sq hj ∇qDðqÞsq ji :(5) Equipped with these electronic and phononic quantities, we move on to the calculation of the various interactions in the electron–phonon system. We start with the e–ph interactions. The vertices (matrix elements) in Bloch space are given by 30 gsmn kq ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ 2mτ_ωsq s1 NeNph X ReRph expðikReþiqRphÞX s0m0n0 gs0m0n0 ReRph Unn0k0Uy m0mkus0sqþgFr;smn kq : (6) where gsmn ReRph are the e–ph matrix elements in Wannier space. The additional final term is the long-range correction due to the socalled Fröhlich interaction 33,34 required for polar materials. For these materials, the gReRph is constructed to be the short-range part of the full interaction. The Fröhlich term takes the following form: gFr;smn kq ¼ie2 ε0VP τffiffiffiffiffiffiffiffiffiffiffiffiffi _ 2mτ_ωsq qUk0Uy k hi nm P G≠q ðqþGÞZ τξτ;sq ðqþGÞϵ1ðqþGÞ ´exp iðqþGÞrτ ½: (7) where ξ τ,sq =u τ,sq is the eigendisplacement of the basis atom τdue to the phonon mode sq,r τ is the position of the atom τin the primitive unit cell, and the sum over the reciprocal lattice vectors G is performed using the Ewald sum technique by multiplying each term of the summation by the factor exp½ðqþGÞϵ1 ðqþGÞ=ð4αÞ, where the parameter αis set to 1 Bohr −2 . Recently, it has been shown that quadrupolar corrections may be important for an accurate description of the e-acoustic phonon interactions 35,36 . Particularly, the electron mobility can be strongly affected for those piezoelectric materials in which the transport-active band extremum is at the BZ center. Currently, we do not have the capabilities for including the quadrupolar corrections. As such, care must be taken when dealing with strongly polar materials such as cubic GaAs and wurtzite GaN that feature zone-centered band extrema. We plan to include support for quadrupolar corrections in a future release of the code. In terms of the e–ph matrix elements described above, the temperature-dependent transition rates of the electrons due to the phonon absorption (+) and emission (−)processesare given by 37 Xe-ph;þ mhkink0jsq Xe-ph; mhkink0jsq 8 < :9 = ;¼2π _Nk gsmn hkiq  2f0 mhkið1f0 nk0Þn0 sqδðϵnk0ϵmhki_ωsqÞ f0 mhkið1f0 nk0Þð1þn0 sqÞδðϵnk0ϵmhkiþ_ωsqÞ () ; (8) where f0(n 0 ) is the equilibrium, i.e., the Fermi (Bose), distribution of the electron (phonon) gas, and the delta functions enforce the energy conservation in a scattering process. The notation hkik0jq above means that the initial electron wave vector is taken to be in the irreducible Brillouin zone (IBZ), the final electron wave vector is on the full first Brillouin zone (FBZ), the phonon wave vector mediating this transition is q¼½k0hki,withthe square brackets denoting a modulo operation with a reciprocal lattice vector. N k is the number of electronic wave vectors in the FBZ. N.H. Protik et al. 2 npj Computational Materials (2022) 28 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences 1234567890():,; In elphbolt, we have the option of including electroncharged impurity (e-chimp) scattering to capture the effect of charged dopants on the transport properties. The e-chimp interaction is calculated within the first Born approximation, taking the impurity potential to have a static, Yukawa (screened Coulomb) form. The modulus squared of the interaction vertex is given by (generalized from ref. 38 ) ge-chimp q  2¼1 VX i ni Zie ε0ϵ0ðq2þq2 TFÞ  2 ;(9) where i=p,ndenotes p-orn-type doping, n i is the concentration of the dopant, Z i is the ionization of the dopant, eis the absolute electronic charge, ϵ 0 is the zero frequency dielectric constant of the material, qis the wave vector magnitude measured from the nearest BZ center, and q TF is the Thomas–Fermi screening wave vector given by q2 TF ¼dse2β NkVε0ϵ0X mk f0 mk1f0 mk  ;(10) where d s is the spin degeneracy of the electronic state, and β ðkBTÞ1in terms of the Boltzmann constant k B and temperature T. In terms of the matrix elements, the charged impurity mediated m〈k〉to nk0transition rates are given by Xe-chimp mhkink0jq¼2π _Nk ge-chimp q  2 1hkik0 kk0  f0 mhki1f0 mhki  δðϵnk0ϵmhkiÞ;(11) where we have multiplied the squared matrix elements with a transport factor corresponding to the in-scattering for elastic processes. Note that this term suppresses forward scattering. At the present, we do not include e–e interactions. Typically, the e–e scattering does not strongly affect the transport in metals 37 . However, in degenerate semiconductors, its effects might be strong 39 . A rigorous treatment of the e–e interaction necessarily requires the calculation of the dynamical screening of the electron gas. Typically, this is done in the random phase approximation which introduces an additional Bosonic system—the plasmons. We wish to include e–e and the related e-plasmon scattering in a future release. Next, we look at the interactions of the phonons. We start with the ph–e interaction. The lowest-order process involves two electrons, and the transition probabilities are given by 37 Yph-e shqijmknk0¼2π _Nq gsmn khqi  2 f0 mkð1f0 nk0Þn0 shqiδðϵnk0ϵmk_ωshqiÞ; (12) where N q is the number of phonon wave vectors in the FBZ and k0¼½kþhqi. This describes the same +process given in Eq. (8). The lowest-order ph–ph interaction vertices describing the q00 !q±q0processes are given by 4 V±;ss0s00 hqiq0q00 ¼X 0 iX jk X αβγ Ψαβγ ijk eα i;shqieβ j;s0±q0eγ k;s00q00 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mimjmk p;(13) where i,j,klabel the atoms in the supercell, the primed sum indicates restriction to the central primitive unit cell, the Greek letters are Cartesian directions, ej;s0q0¼expðirjq0Þξj;s0q0is the eigendisplacement of atom jin the supercell due to the phonon mode s0q0, and Ψ ijk are the third-order anharmonic force constants. We will require only the modulus squared of the above expression. As such, we need only calculate the V − processes, since Vþ;ss0s00 qq0q00  2¼V;ss0s00 qq0q00  2 . In terms of the ph–ph vertices, the anharmonic phonon transition rates can be calculated as 4 W± shqis0q0js00q00 ¼π_ 4Nq V±;ss0s00 hqiq0q00  2 ωshqiωs0q0ωs00q00 ´ðn0 shqiþ1Þn0 sq0þ1 2±1 2  n0 sq00δðωshqi±ωsq0ωsq00Þ: (14) We note that four-phonon interactions have been shown to be important for phonon transport in strongly anharmonic materials and those that feature anomalously weak three-phonon scattering rates for modes that dominate the transport at high temperatures; see ref. 40 and the references therein. We do not currently have the functionality for calculating four-phonon interactions but plan to include it in a future release. In elphbolt, we may also consider lowest-order ph–isotope and ph–substitution defect interactions within the first Born approximation and assuming the defect to be an on-site perturbation and low in concentration. This approximation can be described as the Tamura model 41 . The transition rates for these two-phonon processes are given by Wph-x shqis0q0¼πω2 shqi 2n0 shqið1þn0 sqÞX τ gτξ τ;shqiξτ;s0q0  2 δðωs0q0ωshqiÞ; (15) where x={iso, subs} for isotope and substitution defect scattering, respectively, and gτ¼X t ftτ1Mtτ hMiτ  2 (16) is the mass variance parameter of the host atom τ,tdenotes the type—isotope or substitution—of the guest atom, f tτ is the ratio of type-tguest to type-τhost, and 〈M〉 τ is the average on-site mass. It is worth noting that the above description of the ph-defect interaction is simplistic and can fail for defects which cause extended bond perturbations. There exist advanced diagrammatic techniques that can handle such cases better 42–44 . Such methods are planned for inclusion in a future release of elphbolt. Coupled and decoupled BTEs In elphbolt we may consider the electric (E) and temperature gradient (∇T)fields as the drivers of the electron and the phonon currents. These applied fields cause the distribution functions of the electrons and the phonons, f mk and n sq , respectively, to deviate from their equilibrium forms. In the linear response regime, these are given by fmkf0 mk½1þð1f0 mkÞΨmk nsqn0 sq½1þð1þn0 sqÞΦsq;(17) The deviation functions above of the electrons and the phonons, respectively, can be written as Ψmk¼β∇TImkβEJmk Φsq¼β∇TFsqβEGsq;(18) where I mk (F sq ) is the electron (phonon) response function that measures the deviation from an equilibrium of the occupation of the electron (phonon) state mk(sq) due to the applied ∇Tfield. Similarly, J mk (G sq ) is the electron (phonon) response to the Efield. The coupled electron and phonon BTEs have been given in numerous references, e.g., refs. 21,23–26,37 . Here, we write them in N.H. Protik et al. 3 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences npj Computational Materials (2022) 28 the following form in terms of the response functions: ð19Þ In the equations above, the terms with the superscript 0 are the relaxation time approximation (RTA) terms. These involve a direct coupling to the applied field, which is why in the phonon equation to the Efield, the RTA term is identically zero. Furthermore, the RTA only involves out-scattering processes. Truncating the BTEs up to the RTA term is equivalent to ignoring the in-scattering corrections and the drag effect. Next, the in-scattering corrections are given by the terms with the superscript S. These terms are functionals of the response function of the same species and, as such, they are called the self terms. The inclusion of these terms renders the BTEs (barring the one for G) into a set of decoupled equations which can be solved iteratively. Truncating the BTEs up to the self term, however, still ignores the drag effect. This is because, at this level of the approximation, the electron equations implicitly take phonons to remain in equilibrium, while the phonon equations take electrons to remain in equilibrium. Finally, the terms with the superscript D represent the drag terms. These arefunctionalsoftheresponsefunctionoftheotherspecies. Theinclusionofthesetermsallows an accurate description of the mutual drag effect in the interacting electron–phonon system. At this level, each pair of BTEs must be solved selfconsistently. Furthermore, there exists a thermodynamic restriction known as the Kelvin–Onsager relationship 17 that reflects a deeper coupling of the two pairs of equations for the two fields. In elphbolt, we have devised a fast, iterative scheme that allows us to obtain the full solution of the Eqs. (19), while respecting the thermodynamic restrictions mandated by the Kelvin–Onsager relationship. This is discussed further in the Transport coefficients subsection. Now we provide the expressions for each term in the BTEs in Eq. (19). First, we give the expressions for the electronic equations. The field coupling, RTA terms are given by I0 mk J0 mk  ¼ðϵmkμcÞ=T e  vmk WRTA mk ;(20) where μ c is the chemical potential of the electronic system. The electron RTA scattering rates are given by WRTA mk¼1 f0 mkð1f0 mkÞ"P nsq Xþ mknk0jsqþX mknk0jsq  þP nq Xe-chimp mknk0jq#: (21) The summation above of the independent scattering channels at the RTA level is known as Matthiessen’s Rule. The self terms are given by ΔIS mk ΔJS mk  ¼1 f0 mkð1f0 mkÞ 1 WRTA mkX snk0 Ink0 Jnk0  Xþ mknk0jsqþX mknk0jsq  : (22) Lastly, the drag terms are ΔID mk ΔJD mk  ¼1 f0 mkð1f0 mkÞ 1 WRTA mkX snk0 FsqXþ mknk0jsqþFsqX mknk0jsq GsqXþ mknk0jsqþGsqX mknk0jsq 8 < :9 = ; : (23) Similarly, for the phonons we have F0 sq G0 sq () ¼_ωsq=T 0  vsq WRTA sq ;(24) where the phonon RTA scattering rates are given by WRTA sq¼1 n0 sqð1þn0 sqÞ"P s0q0s00q00 Wþ sqs0q0js00q0þ1 2W sqs0q0js00q0  þdsP mnk Ysqjmknk0þP s0q0 Wphx sqs0q0#: (25) The self and drag terms are, respectively ΔFS sq ΔGS sq  ¼1 n0 sqð1þn0 sqÞ 1 WRTA sq ´P s0q0s00q00 "Wþ sqs0q0js00q00 Fs00q00Fs0q0 Gs00q00Gs0q0 no þ1 2W sqs0q0js00q00 Fs00q00þFs0q0 Gs00q00þGs0q0 no #(26) and ΔFD sq ΔGD sq () ¼ds n0 sqð1þn0 sqÞ 1 WRTA sqX mnk Ysqjmknk0 Ink0Imk Jnk0Jmk  :(27) With these expressions, we solve the coupled e–ph BTEs, Eqs. (19), using an iterative procedure which is an improvement over the one used earlier in refs. 25,26 . Here we use a unified scheme, indexed by a single iterator, that of the ph BTE. In this approach, for a coupled BTEs solution, the e BTE is internally iterated to selfconsistency for each iteration of the ph BTE. Thus, once the ph BTE has achieved self-consistency, the e BTE is guaranteed to do so also. This procedure is physically justified since the electron system is, in general, faster than the phonon system. Moreover, this scheme allows the strict enforcement of the Kelvin–Onsager relationship after each ph BTE iteration because the electron system is always brought to consistency with the current nonequilibrium phonon system. However, since the Kelvin–Onsager relationship mandates a thermodynamic constraint over the thermoelectric transport coefficients of both the electron and the phonon systems (see subsection Transport coefficients), in order to confirm that this relationship is satisfied, the coupled BTEs iteration must be performed for both the Eand ∇Tfields, simultaneously. Thus, all four BTEs in Eqs. (19) are tightly coupled. This approach is an improvement over the simple, dual iterator scheme used in refs. 25,26 , both conceptually and numerically. We offer four different levels of solutions of the BTE for each species and field: RTA, partially decoupled, dragless full, and dragged full. Note that the RTA and dragless full solutions for the case of the phonons under the influence of the electric field are trivially zero. In our iterative approach, the computational time and space difference between the partially decoupled solution and the dragged full solution is not large, and both the RTA (iteration 0) and partially decoupled (iteration 1) solutions are provided en route to the dragged full solution. By default, elphbolt will carry out a set of decoupled, hence, dragless, ph and e BTEs immediately after the coupled BTEs are fully solved. This allows the users to compare the transport coefficients in all the abovementioned approximation levels after a single computational run. The users, however, also have the option of calculating N.H. Protik et al. 4 npj Computational Materials (2022) 28 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences only the decoupled e and ph BTEs (dragless full), if they wish to do so. Transport coefficients From the solutions of the response functions, we can obtain the following transport coefficient tensors. From the electronic charge current, we get σ σS no ¼dse VkBTX mk f0 mkð1f0 mkÞvmkJmk Imk  ;(28) where σis the electronic charge conductivity and Sis the (Seebeck) thermopower. From the electronic heat current, we get αel κ0;el  ¼ ds VkBTX mkðϵmkμcÞf0 mkð1f0 mkÞvmkJmk Imk  ; (29) where α el is closely related to the electronic Peltier coefficient. The electronic component of the (Peltier) thermopower is given by Q el =α el (σT) −1 . The tensor κ 0,el is the electronic thermal conductivity in the zero Efield (closed circuit) condition. The electronic thermal conductivity in the open circuit condition, which is what can be measured in experiments, is given by κ el =κ 0,el −α el S. Lastly, from the phonon heat current, we obtain αph κph  ¼1 VkBTX sq _ωsqn0 sqð1þn0 sqÞvsqGsq Fsq  ;(30) where κ ph is the phonon thermal conductivity and α ph is related to the phonon Peltier coefficient. Specifically, the phonon component of the (Peltier) thermopower is given by Q ph =α ph (σT) −1 . The Kelvin–Onsager relationship unifies the Seebeck and the Peltier effects and mandates that σS¼αT1;(31) where α=α el +α ph17 . Note that this implies that the same thermopower Q≡S=Q el +Q ph is obtained by both the Seebeck and the Peltier pictures. As such, the above relationship connects the ∇Tfield e BTE and the ∇Tand Efield ph BTEs. Since, for a given field, the BTE foronespeciesiscoupledtotheonefortheotherviathedrag term, the Kelvin–Onsager relationship effectively leads to a tight physical connection between all the four BTEs of the interacting e–ph system. However, during the course of the iterations, numerical issues may cause the system of equations to deviate from the Kelvin–Onsager relationship. To remedy this, corrective measures must be employed. Now, within the iterative scheme described in the previous subsection, it is straightforward to enforce the Kelvin–Onsager relationship. To do this, we split σSinto a diffusion and a drag term and draw a parallel with the Peltier decomposition: σS¼σSdiff þσSdrag αT1¼αelT1þαphT1:(32) Then, decomposing Ias (band and wave vector indices dropped for brevity) I¼Idiff þIdrag (33) and demanding σS diff [I diff ]=α el T −1 , we can arrive at the relation: Idiff ¼ϵμc eT J:(34) The problem of enforcing the Kelvin–Onsager relationship then boils down to satisfying σSdrag½λIdrag¼αphT1;(35) where λis a small, corrective scalar which can be easily found using a bisection method. Implementation The coupled BTEs solver described above is implemented in Fortran 2018. This allows us to make use of the object-oriented programming (OOP) support and the built-in coarray functionality that provides concise, the native syntax for parallelization. Specifically, we create the following seven derived types: crystal,symmetry,numerics,electron,phonon,epw_wannier, and bte dealing with the components of the problem that the names suggest. Each derived type contains its own data and procedures (functions and subroutines). Apart from these, there are separate modules for immutable parameters, helper procedures, etc. This hybrid OOP/procedural design enables the extensibility of the code. Boilerplate getter and setter functions are generally avoided. Instead, the intent and use, only keywords of Fortran are strictly used to control the read, write, and use access of the different components of the code. This design strategy makes the code compact (about 6700 lines) for what it offers and easily readable. As a general rule, code repetition is avoided unless the generalizations lead to a slow or physically unclear source. We tried to strike a balance between the speed of development, execution, readability, and extensibility. Below we discuss the general workflow and the logical structure of the program. Workflow and structure In Fig. 1we present the workflow of elphbolt.Forafull calculation complete with the e–ph drag effect, we require the second-order interatomic force constants (IFC2s) from Quantum Espresso 45–47 ; the third-order interatomic force constants (IFC3s) from thirdorder.py 4 which provides an interface with Quantum Espresso; and the Wannier space electronic Hamiltonian, dynamical matrix, e–ph matrix elements, and the real space cell maps and degeneracies from EPW 5,30,34 .EPW internally interfaces with Quantum Espresso and Wannier90 48 . The generation of the IFC2s and IFC3s are also part of the ShengBTE workflow, and the users of that code will find that only one extra step—i.e., an EPW calculation—is required to generate all the input data for an e–ph coupled BTEs calculation with elphbolt. We have provided two modified EPW source files that allow the generation of some of the required Wannier space information. The user also needs to provide an input file called input.nml.Theformatofthe input is described in details in the file README.org. Next, we outline the logical flow of the program. The main program is in the file elphbolt.f90. The program starts with Fig. 1 Workflow of elphbolt showing the various input data required for a coupled e–ph BTEs calculation. Data generated by the external codes shown in the boxes on the left are passed into elphbolt, which carries out the transport calculations to generate the results listed on the right. N.H. Protik et al. 5 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences npj Computational Materials (2022) 28 creating objects of the seven derived types mentioned earlier. Then, following a welcome message, it initializes the crystal object. This involves reading the information about the crystal and initializing the appropriate internal variables. Following this, the numerics object is initialized which involves the reading in of the transport wave vector meshes, data output directory, BTE solution type, etc. Next, the symmetry object is initialized. At this step, the symmetries of the crystal and the BZ are calculated. Following this, the epw_wannier object is initialized which involves reading in the Wannier space information. Next, we initialize the electron object. In this step, the information about the bands, transport-active energy window, chemical potential, etc. are read in. The transport energy window restricted electronic IBZ is generated along with the bands, eigenstates, and velocities. Next, we initialize the phonon object which involves the calculation of the phonon branches, eigenstates, and velocities. The next major step is the calculation of the interaction vertices. These zero temperature quantities are stored in the disk for later temperature and carrier concentration sweeps. The temperaturedependent transition rates are calculated and stored in the disk also for their reuse in the various types of BTE solutions that are available. The transition rates expressions include the energyconserving delta functions. We provide two methods for the evaluation of these delta functions: the triangular method 49,50 and the analytic tetrahedron method 51 . The first is the only option for 2d systems and is the default for 3d systems. Unlike the more commonly used Gaussian or Lorentzian methods, the triangular and tetrahedron methods do not have any additional smearing parameter, and, as such, the electron and phonon wave vector meshes are the only parameters to converge. Finally, the bte object is used to solve the BTEs. The code produces a large amount of data for analysis, including the RTA scattering rates, band/branch resolved transport tensors, response functions at the RTA, partially decoupled, and fully iterated levels, among others. The users may also run the code in a post-processing mode to generate spectral transport coefficients if required. Full descriptions of all the input options and the output data are provided in the README.org file. Example: cubic silicon As a demonstration of the code, we calculate the drag effect on the thermopower in cubic silicon. In addition, we include calculations of the mobility and thermal conductivity. We consider first an n-type doped sample of low carrier density 2.75 × 10 14 cm −3 over a range of temperatures. This is close to the carrier concentration in the best sample (number 537) of Geballe and Hull’s experimental work 15 . We show in Fig. 2a comparison of the calculated and the measured thermopower magnitudes. The experimental data is collected from Fig. 1 of ref. 15 . The black curve is the calculated total thermopower. In the Seebeck picture, this includes the phonon drag contribution. And in the Peltier picture, this is the sum of the electronic contribution and the phonon contribution, the latter being purely due to the electron drag effect. The Peltier picture also allows a clean separation of the electronic and the phonon contributions. This Peltier breakdown of the total thermopower is shown in Fig. 2. The solid blue curve is the electronic contribution to the thermopower. This contribution decreases slightly with decreasing temperature. The solid green curve is the phonon thermopower. This starts off as significantly lower than the electronic thermopower at 300 K, but quickly overtakes the latter around 175 K, before dominating by more than an order of magnitude at 50 K. The calculated total thermopower is in excellent agreement with the experimental measurements (red circles). Without the drag effect, the phonon contribution would be trivially zero, and the electronic contribution would be an order of magnitude off from the experimental values at low temperatures. Furthermore, without drag, the temperature dependence of the thermopower would be spectacularly wrong. The explanation of the origin of the strong drag effect in silicon has previously been given in refs. 21,23 , and is not reproduced here. These earlier ab initio works captured the strong drag behavior using a partially decoupled solution of the e and ph BTEs. Nevertheless, they found excellent agreement with experimental measurements. Our calculations corroborate their finding that the partially decoupled solution is indeed sufficient to capture the strong drag effect in silicon—the green and blue squares denote the partially decoupled solutions, and they coincide nearly perfectly with corresponding full drag curves. For this case, the RTA and dragless full solutions of the e BTE also give essentially the same result and are not shown on the plot to reduce clutter. However, this is not guaranteed to hold true for all materials, and, in general, a full drag solution of the e–ph BTEs should be used. Figure 3shows the calculated thermopower of an n-type sample with a high density of 2.7 × 10 19 cm −3 , matching that of sample 140 in ref. 15 . Good agreement with measured data is again obtained. The percentage difference of the calculated total thermopower from the experimental value near 50 (300) K is around 15 (1)%. Note that the calculated thermopower without Fig. 2 Temperature dependence of the thermopower of silicon for an n-type carrier concentration of 2.75 × 10 14 cm −3 .The red circles are measurements (sample 537, concentration 2.8 × 10 14 cm −3 )by Geballe and Hull 15 . Fig. 3 Temperature dependence of the thermopower of silicon for an n-type carrier concentration of 2.7 × 10 19 cm −3 .The red circles are measurements (sample 140, concentration 2.7 × 10 19 cm −3 )by Geballe and Hull 15 . N.H. Protik et al. 6 npj Computational Materials (2022) 28 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences the phonon component again significantly underestimates the measured thermopower across the full range of temperatures considered. In Fig. 4, we present the calculated temperature-dependent mobilities for a carrier concentration of 2.75 × 10 14 cm −3 . Excellent agreement is found over the entire temperature range considered here with experimental measurements (shown in red circles) on similar low-doped samples 52 . The effect of phonon drag on this quantity is negligible. In fact, the difference between the dragged full, RTA, dragless full, and partially decoupled solutions is very small, and, to reduce clutter, we do not show the latter two results on the plot. This corroborates the finding in ref. 21 . Figure 5shows the calculated temperature-dependent mobilities for a degenerate carrier concentration of 2 × 10 19 cm −3 and the comparison to experimental data from samples with similar electron concentrations from ref. 53 . The calculated mobilities including e-chimp scattering (solid blue curve) are between a factor of 4 to 5 higher than the measured values (red circles). The RTA, dragless full, and partially decoupled solutions give nearly the same values as the dragged full solution and, to reduce clutter, we do not show these points on this plot. The large discrepancy between the calculated and the measured mobilities could be due to a combination of multiple reasons. First, the measurements were done on compensated samples and the acceptor and donor compositions were not reported in ref. 53 . Our calculations have assumed that the mobile electron concentration is equal to the ionized donor density, while the acceptor density has been taken to be zero. In the actual samples, acceptors provide additional charge scattering centers, thus lowering the mobility. Second, the e-chimp scattering employed in the calculation assumes a static, Thomas–Fermi screened Coulomb interaction treated in the Born approximation. It is known that the Thomas–Fermi model leads to a significant overscreening of the interaction in the degenerate limit 38 . The validity of the Born approximation is also stretched in this limit and more sophisticated non-perturbative approaches might be better suited. Lastly, the consideration of the e-plasmon scattering, which is currently not included in our calculation, has been shown to reduce the mobility of silicon significantly at high carrier concentrations 39 .A more rigorous treatment of the e-chimp scattering along with e–e scattering is planned for a future version of elphbolt. We also show in Fig. 5the calculated mobility with the e-chimp interaction turned off (dashed blue curve). Here we envision that the mobile charge carriers are created in a region in which charged dopants do not exist, which could happen, for example, through modulation doping 54 , or in GaN/AlGaN heterostructures 55 . In this case, the phonon drag effect leads to a large increase in mobility—at 300 (50) K, the drag gain of mobility is about a factor of 2.5 (50). Such high gains in mobility can, potentially, be exploited along with those previously identified in the thermopower 21,23,56 to boost the thermoelectric figure-ofmerit. Note that, in this case, the phonon drag effect is not fully captured by the partially decoupled solution, which predicts about half the value given by the dragged full solution at 50 K. At this temperature, the RTA and dragless full solutions undershoot the value of the dragged full solution by about a factor of 5. Finally, in Fig. 6we compare the calculated phonon thermal conductivities against measurements on high purity samples of natural silicon. The measured values (red circles) are taken from ref. 57 . Remarkable agreement is found over the full temperature range considered for the low concentration case (solid blue curve). We also plot the results for a high concentration case in dashed green. The ph–e interactions cause a weaker suppression at high temperatures compared to that in the low-temperature limit even at a high carrier concentration of 2 × 10 19 cm −3 . This happens Fig. 4 Temperature dependence of the mobility of silicon for an n-type carrier concentration of 2.75 × 10 14 cm −3 .The red circles are measurements on various different samples with carrier concentrations ranging from 3.5 × 10 13 to 1.4 × 10 14 cm −352 . Fig. 5 Temperature dependence of the mobility of silicon for an n-type carrier concentration of 2 × 10 19 cm −3 .The red circles are measurements by Yamanouchi et. al. 53 on samples 23, 25, and 32, with dopant densities 2.28 × 10 19 , 2.42 × 10 19 , and 3.35 × 10 19 cm −3 , respectively. Fig. 6 Temperature dependence of the phonon thermal conductivity of silicon for n-type carrier concentrations of 2.75 × 10 14 and 2 × 10 19 cm −3 .The red circles are measurements on high purity samples with natural isotopic mix 57 . N.H. Protik et al. 7 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences npj Computational Materials (2022) 28 because the ph–ph scattering rates dominate over the ph–e ones at high temperatures. This is consistent with the findings in ref. 58 . At low temperatues, the low energy acoustic phonons progressively contribute more to the thermal conductivity, while, at the same time, the ph–e scattering rates begin to dominate over the ph–ph ones for these modes. This results in a stronger suppression of the thermal conductivity at low temperatures for the high doped case. The electron drag effect on the phonon thermal conductivity is found to be small in this material, as has been shown earlier in refs. 21,23 . This has also been shown to be true for gallium arsenide 25 and silicon carbide 26 . The reason for this has been discussed in these references. The dragged full, dragless full, partially decoupled, and the RTA solutions all give nearly the same results and only the first type is shown on the plot. DISCUSSION In this work, we discussed the theory and implementation of elphbolt—an efficient code for solving the coupled electron–phonon Boltzmann transport equations. This gives ab initio access to the thermal, charge, and thermoelectric transport properties in materials, and the effect of e–ph drag on them. The code is distributed as Free/Libre software under the GNU General Public License version 3 that gives the user the right to use, modify, and distribute the original and their modified versions of the software. The code combines concepts of the object-oriented and procedural programming styles, has a clean, modular structure, features coarray parallelization, and is well-documented. This makes the code easy to extend. For future releases, we plan to include a more rigorous treatment of impurity scattering, electron–electron interactions, provisions for including four-phonon interactions 59 ,quadrupolar corrections, and magnetotransport. METHODS Calculation details We use the norm-conserving, Perdew–Zunger, local density approximation (LDA) pseudopotential 60 . A relaxed lattice constant of 5.40 Åis found. The phonon calculation is performed using a 12 × 12 × 12 k-mesh and a 6 × 6 × 6q-mesh. The EPW calculation is done with four valence and four conduction bands, using an initial guess of sp 3 projections on the two basis atoms. The IFC3s are calculated using a 5 × 5 × 5 supercells (250 atoms) with a six nearest neighbor cut-off and Γ-point sampling. In the transport calculations, we include e–ph, e-chimp (assuming singly charged impurities), ph–ph, ph–e, and ph-iso scattering. For the transport calculations, converged 50 × 50 × 50 qand 150 × 150 × 150 k-meshes, denoted (50, 150), are used. We use a gcc build of elphbolt with the OpenCoarrays library 61 for coarray support. A typical calculation of the fully coupled e–ph BTEs using the (50, 150) mesh set takes on the order of 3000 CPU-hours. A breakdown of the various important components is given in Table 1. These numbers were calculated on four nodes each equipped with 28 Intel(R) Xeon(R) CPU E52680 v4 @ 2.40 GHz cores for an n-type carrier concentration of 2.75 × 10 14 cm −3 at 300 K and a relative convergence threshold of 0.0001. The coupled e–ph BTEs required six iterations to converge, whereas the decoupled e and ph BTEs took four and seven iterations, respectively. The total run time, of course, can vary significantly depending on the speed of disk read/write and that of the CPUs, the temperature, and the carrier concentration. Note that once the interaction vertices are calculated, these can then be reused during the solutions of the BTEs for different temperatures and carrier concentrations. DATA AVAILABILITY The input files needed to generate both the force constants and Wannier space data required to reproduce the results in this work are available from github 27 . CODE AVAILABILITY The code used in this work is available from github 27 . Received: 24 September 2021; Accepted: 11 January 2022 Published online: 07 February 2022 REFERENCES 1. Hohenberg, P. & Kohn, W. Inhomogeneous electron gas. Phys. Rev. 136, B864 (1964). 2. Kohn, W. & Sham, L. J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133 (1965). 3. Baroni, S., De Gironcoli, S., Dal Corso, A. & Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 73, 515 (2001). 4. Li, W., Carrete, J., Katcho, N. A. & Mingo, N. ShengBTE: a solver of the Boltzmann transport equation for phonons. Comput. Phys. Commun. 185, 1747–1758 (2014). 5. Poncé, S., Margine, E. R., Verdi, C. & Giustino, F. EPW: electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions. Comput. Phys. Commun. 209, 116–133 (2016). 6. Zhou, J.-J. et al. Perturbo: a software package for ab initio electron–phonon interactions, charge transport and ultrafast dynamics. Comput. Phys. Commun. 264, 107970 (2021). 7. Broido, D. A., Malorny, M., Birner, G., Mingo, N. & Stewart, D. A. Intrinsic lattice thermal conductivity of semiconductors from first principles. Appl. Phys. Lett. 91, 231922 (2007). 8. Liu, T.-H., Zhou, J., Liao, B., Singh, D. J. & Chen, G. First-principles mode-by-mode analysis for electron-phonon scattering channels and mean free path spectra in GaAs. Phys. Rev. B 95, 75206 (2017). 9. Poncé, S., Margine, E. R. & Giustino, F. Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors. Phys. Rev. B 97, 121201 (2018). 10. Bloch, F. Zum elektrischen Widerstandsgesetz bei tiefen Temperaturen. Z. f.ür. Phys. 59, 208–214 (1930). 11. Peierls, R. Zur Theorie der elektrischen und thermischen Leitfähigkeit von Metallen. Ann. Phys. 396, 121–148 (1930). 12. Gurevich, Y. G. & Mashkevich, O. L. The electron-phonon drag and transport phenomena in semiconductors. Phys. Rep. 181, 327–394 (1989). 13. Frederikse, H. P. R. Thermoelectric power of germanium below room temperature. Phys. Rev. 92, 248 (1953). 14. Geballe, T. H. & Hull, G. W. Seebeck effect in germanium. Phys. Rev. 94, 1134 (1954). 15. Geballe, T. H. & Hull, G. W. Seebeck effect in silicon. Phys. Rev. 98, 940 (1955). 16. Herring, C. Theory of the thermoelectric power of semiconductors. Phys. Rev. 96, 1163 (1954). 17. Sondheimer, E. H. The Kelvin relations in thermo-electricity. Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 234, 391–398 (1956). 18. Cantrell, D. G. & Butcher, P. N. A calculation of the phonon-drag contribution to the thermopower of quasi-2D electrons coupled to 3D phonons. I. General theory. J. Phys. C. Solid State Phys. 20, 1985 (1987). 19. Cantrell, D. G. & Butcher, P. N. A calculation of the phonon-drag contribution to the thermopower of quasi-2D electrons coupled to 3D phonons. II. Applications. J. Phys. C. Solid State Phys. 20, 1993 (1987). 20. Mahan, G. D., Lindsay, L. & Broido, D. A. The Seebeck coefficient and phonon drag in silicon. J. Appl. Phys. 116, 245102 (2014). 21. Zhou, J. et al. Ab initio optimization of phonon drag effect for lower-temperature thermoelectric energy conversion. Proc. Natl Acad. Sci. USA 112, 14777–14782 (2015). 22. Zhou, J. et al. First-principles simulation of electron transport and thermoelectric property of materials, including electron-phonon scattering, defect scattering, and phonon drag. Mater. Cloud Arch.2020 (2020). Table 1. Breakdown of the computational time needed for the various expensive parts of the code for silicon. An n-type carrier concentration of 2.75 × 10 14 cm −3 at 300 K using a (50, 150) mesh set is considered here. task e–ph vertex ph–ph vertex e–ph BTEs e BTE ph BTE time (CPU-hours) 744 545 1467 66 28 N.H. Protik et al. 8 npj Computational Materials (2022) 28 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences 23. Fiorentini, M. & Bonini, N. Thermoelectric coefficients of n-doped silicon from first principles via the solution of the Boltzmann transport equation. Phys. Rev. B 94, 85204 (2016). 24. Macheda, F. & Bonini, N. Magnetotransport phenomena in p-doped diamond from first principles. Phys. Rev. B 98, 201201 (2018). 25. Protik, N. H. & Broido, D. A. Coupled transport of phonons and carriers in semiconductors: a case study of n-doped GaAs. Phys. Rev. B 101, 75202 (2020). 26. Protik, N. H. & Kozinsky, B. Electron-phonon drag enhancement of transport properties from a fully coupled ab initio Boltzmann formalism. Phys. Rev. B 102, 245202 (2020). 27. Protik, N. H. elphbolt. https://github.com/nakib/elphbolt (2021). 28. Marzari, N. & Vanderbilt, D. Maximally localized generalized Wannier functions for composite energy bands. Phys. Rev. B 56, 12847 (1997). 29. Souza, I., Marzari, N. & Vanderbilt, D. Maximally localized Wannier functions for entangled energy bands. Phys. Rev. B 65, 35109 (2001). 30. Giustino, F., Cohen, M. L. & Louie, S. G. Electron-phonon interaction using Wannier functions. Phys. Rev. B 76, 165108 (2007). 31. Feynman, R. P. Forces in molecules. Phys. Rev. 56, 340 (1939). 32. Pick, R. M., Cohen, M. H. & Martin, R. M. Microscopic theory of force constants in the adiabatic approximation. Phys. Rev. B 1, 910 (1970). 33. Sjakste, J., Vast, N., Calandra, M. & Mauri, F. Wannier interpolation of the electronphonon matrix elements in polar semiconductors: polar-optical coupling in GaAs. Phys. Rev. B 92, 54307 (2015). 34. Verdi, C. & Giustino, F. Fröhlich electron-phonon vertex from first principles. Phys. Rev. Lett. 115, 176401 (2015). 35. Brunin, G. et al. Electron-phonon beyond Fröhlich: dynamical quadrupoles in polar and covalent solids. Phys. Rev. Lett. 125, 136601 (2020). 36. Jhalani, V. A., Zhou, J.-J., Park, J., Dreyer, C. E. & Bernardi, M. Piezoelectric electronphonon interaction from ab initio dynamical quadrupoles: impact on charge transport in wurtzite GaN. Phys. Rev. Lett. 125, 136602 (2020). 37. Smith, H. & Jensen, H. H. Transport Phenomena (Clarendon Press, Oxford Univ. Press, 1989). 38. Chattopadhyay, D. & Queisser, H. J. Electron scattering by ionized impurities in semiconductors. Rev. Mod. Phys. 53, 745 (1981). 39. Caruso, F. & Giustino, F. Theory of electron-plasmon coupling in semiconductors. Phys. Rev. B 94, 115208 (2016). 40. Feng, T. & Ruan, X. In Nanoscale Energy Transport (ed. Liao, B.) Ch.2 (IOP, 2020). 41. Tamura, S.-i Isotope scattering of dispersive phonons in Ge. Phys. Rev. B 27, 858 (1983). 42. Mingo, N., Esfarjani, K., Broido, D. A. & Stewart, D. A. Cluster scattering effects on phonon conduction in graphene. Phys. Rev. B 81, 45408 (2010). 43. Katcho, N. A., Carrete, J., Li, W. & Mingo, N. Effect of nitrogen and vacancy defects on the thermal conductivity of diamond: an ab initio Green’s function approach. Phys. Rev. B 90, 94117 (2014). 44. Fava, M. et al. How dopants limit the ultrahigh thermal conductivity of boron arsenide: a first principles study. npj Comput. Mater. 7,1–7 (2021). 45. Giannozzi, P. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys. Condens. Matter 21, 395502 (2009). 46. Giannozzi, P. et al. Advanced capabilities for materials modelling with quantum ESPRESSO. J. Phys. Condens. Matter 29, 465901 (2017). 47. Giannozzi, P. et al. Quantum ESPRESSO toward the exascale. J. Chem. Phys. 152, 154105 (2020). 48. Pizzi, G. et al. Wannier90 as a community code: new features and applications. J. Phys. Condens. Matter 32, 165902 (2020). 49. Kurganskii, S. I., Dubrovskii, O. I. & Domashevskaya, E. P. Integration over the twodimensional Brillouin zone. Phys. Status Solidi 129, 293–299 (1985). 50. Wang, T., Carrete, J., van Roekeghem, A., Mingo, N. & Madsen, G. K. H. Ab initio phonon scattering by dislocations. Phys. Rev. B 95, 245304 (2017). 51. Lambin, P. & Vigneron, J.-P. Computation of crystal Green’s functions in the complex-energy plane with the use of the analytical tetrahedron method. Phys. Rev. B 29, 3430 (1984). 52. Canali, C., Jacoboni, C., Nava, F., Ottaviani, G. & Alberigi-Quaranta, A. Electron drift velocity in silicon. Phys. Rev. B 12, 2265 (1975). 53. Yamanouchi, C., Mizuguchi, K. & Sasaki, W. Electric conduction in phosphorus doped silicon at low temperatures. J. Phys. Soc. Jpn. 22, 859–864 (1967). 54. Störmer, H. L., Dingle, R., Gossard, A. C., Wiegmann, W. & Sturge, M. D. Twodimensional electron gas at a semiconductor-semiconductor interface. Solid State Commun. 29, 705–709 (1979). 55. Ambacher, O. et al. Two-dimensional electron gases induced by spontaneous and piezoelectric polarization charges in N-and Ga-face AlGaN/GaN heterostructures. J. Appl. Phys. 85, 3222–3233 (1999). 56. Yalamarthy, A. S. et al. Significant phonon drag enables high power factor in the AlGaN/GaN two-dimensional electron gas. Nano Lett. 19, 3770–3776 (2019). 57. Inyushkin, A. V., Taldenkov, A. N., Gibin, A. M., Gusev, A. V. & Pohl, H.-J. On the isotope effect in thermal conductivity of silicon. Phys. Status Solidi 1, 2995–2998 (2004). 58. Liao, B. et al. Significant reduction of lattice thermal conductivity by the electronphonon interaction in silicon with high carrier concentrations: a first-principles study. Phys. Rev. Lett. 114, 115901 (2015). 59. Han, Z., Yang, X., Li, W., Feng, T. & Ruan, X. FourPhonon: an extension module to ShengBTE for computing four-phonon scattering rates and thermal conductivity. Comput. Phys. Commun. 270, 108179 (2022). 60. Perdew, J. P. & Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B 23, 5048 (1981). 61. Fanfarillo, A. et al. OpenCoarrays: open-source transport layers supporting coarray Fortran compilers. In Proc. 8th International Conference on Partitioned Global Address Space Programming Models.1–11 (ACM, 2014). ACKNOWLEDGEMENTS This project was funded by the EU-H2020 through H2020-NMBP-TO-IND project GA n. 814487 (INTERSECT). ICN2 is supported by the Severo Ochoa program from Spanish MINECO (Grant No. SEV-2017-0706) and the CERCA Program of Generalitat de Catalunya. Work at Boston College (contributions to code testing and ab initio thermoelectric transport calculations for silicon) was supported by the US Department of Energy (DOE), Office of Science, Basic Energy Sciences under award # DE-SC0021071. NHP acknowledges helpful discussions with Vladimir Dikan, José María Escartín, Xavier Cartoixà, and Riccardo Rurali. We thankfully acknowledge the computer resources at MareNostrum and La Palma and the technical support provided by Barcelona Supercomputing Center (FI-2021-1-0016) and the Center for Astrophysics in La Palma (QS-2021-1-0022), respectively. We also acknowledge computational support from the Boston College Linux clusters and those at ICN2 provided by Grant PGC2018-096955-B-C43 funded by MCIN/AEI/10.13039/ 501100011033 and “ERDF A way of making Europe”. AUTHOR CONTRIBUTIONS NHP is the developer and maintainer of the elphbolt code. NHP performed some of the ab initio calculations and wrote the first draft of the manuscript in consultation with and under the supervision of MP and PO. CL performed some of the ab initio calculations under the supervision of DB. All authors took part in the preparation of the manuscript. COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION Correspondence and requests for materials should be addressed to Nakib H. Protik or Pablo Ordejón. Reprints and permission information is available at http://www.nature.com/ reprints Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/. © The Author(s) 2022, corrected publication 2021 N.H. Protik et al. 9 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences npj Computational Materials (2022) 28