DFTBephy: A DFTB-based approach for electron–phonon coupling calculations

The calculation of the electron–phonon coupling from first principles is computationally very challenging and remains mostly out of reach for systems with a large number of atoms. Semi-empirical methods, like density functional tight binding (DFTB), provide a framework for obtaining quantitative results at moderate computational costs. Herein, we present a new method based on the DFTB approach for computing electron–phonon couplings and relaxation times. It interfaces with phonopy for vibrational modes and dftb+ to calculate transport properties. We derive the electron–phonon coupling within a non-orthogonal tight-binding framework and apply them to graphene as a test case.


Introduction
The phonon-limited mobility is often a critical performance parameter for the application of novel materials in electronic devices.Its computational modeling can be quite challenging since it requires the knowledge of the electronic structure, the phonon modes, and the coupling strength of electrons and phonons [1].Therefore, modeling of the mobility is often based on approximations, such as the effective mass or the deformation potential [2].On the other hand, over the past decade, the ab initio treatment of electron-phonon couplings became much more feasible due to continuous developments of computational methods (e.g., in ABINIT [3], ATK [4], Quantum ESPRESSO [5,6], EPW [7], PER-TURBO [8]) and steadily increasing hardware resources.
For systems with a large number of atoms in the unit cell, such as covalent organic frameworks (COFs) [9], using ab initio approaches is still challenging.Especially in situations where high-throughput screening of many such materials is required, alternative methods are needed.Density functional tight binding (DFTB) [10] is such an approach as it effectively reduces the complexity of density functional theory (DFT), casting the Kohn-Sham equations into a tightbinding form.The method is now rich in extensions [11] and has been successfully applied to study the electronic and structural properties of a vast range of materials.A nonexhaustive list includes organic polymers, COFs [12] and bio-molecular systems [13], transition metal oxides (TiO 2 [14], ZnO [15]), MoS 2 films and nanostructures [16], graphene defects [17] and carbon allotropes.It has been applied particularly to study the structural and electronics of nanoparticles and nanorods of several inorganic materials (Si, SiC, Ag, Au, Fe, Mg) with unfeasible sizes for DFT calculations.Green's function extensions of DFTB have been used to study electron and phonon transport in different nanostructures in the ballistic regime [18].
In this contribution, we focus on the Boltzmann transport theory in the relaxation-time approximation.For this, we first derive an expression for electron-phonon couplings starting from a general non-orthogonal tight-binding Hamiltonian.As such, our result is applicable to DFTB and other tight-binding approaches such as xTB [19].We discuss some subtleties arising in the derivation.The calculation of the electron-phonon couplings as well as the computation of relaxation times and transport properties has been implemented in a Python package dftbephy driving DFTB calculations.As an exemplary demonstration, we consider graphene which has been widely studied in past decades [1,20,21].

The DFTB approach
In DFTB the total energy is expanded in a Taylor series around a suitable reference density [22], 0 , where the individual contributions are given by here E ion−ion and E H denote the electrostatic interaction of ions and electrons (Hartree), respectively.The exchangecorrelation energy and potential are indicated by the index xc and the effective Kohn-Sham single-particle Hamiltonian is expressed as, Starting from this expansion, several DFTB flavors are available.The starting ground is the so-called non-selfconsistent approximation (DFTB1) where only the first two terms in Eq. (1a) are retained.The single-particle wave function, Ψ i (⃗ r) , is expanded in terms of atomic orbitals which are computed as solutions of isolated atoms with a compressed potential.The resulting Hamiltonian matrix elements in Eq. (1c), which are multiplied by the occupation factor f i , neglect three center integrals, crystal field, and pseudo-potential contributions, leaving the evaluation of two center integrals of atomic dimers that can be pre-tabulated.The argument for neglecting three center and crystal field terms is the partial cancelation with the pseudo-potential (2) term ensuring orthogonality of the valence orbitals to all core states, a contribution that vanishes in the two-center approximation [23].The term E 0 [ 0 ] accounts for intera- tomic repulsion and is usually obtained from reference DFT calculations rather than an explicit evaluation of Eq. (1b).
The so-called self-consistent charge approach (SCC-DFTB or DFTB2) includes the second-order term, Eq. (1d), traditionally expanding (⃗ r) in terms of spherical atomic charge densities and Mulliken charges.Coulomb contributions enter both as local terms due to perturbations of the self-consistent Hamiltonian as well as longrange scattering due to macroscopic electric fields, typical of polar crystals.The former is included when performing SCC-DFTB calculations of the electron and phonon band structures and when computing the electron-phonon couplings.The long-range Coulomb scattering is typically added by means of a Fröhlich Hamiltonian, but this is not included in this work.

Electron-phonon couplings
In this section, we give a brief summary of the derivation of the electron-phonon coupling in general non-orthogonal tightbinding Hamiltonians.Several derivations can be found in the literature [20,[24][25][26], but since not all agree with each other in the final result, we think it is useful to discuss our derivation.
We start with explicit reference to a local basis set, � ⟩ , and the dual (biorthogonal) basis, � ⟩ , such that ⟨ � ⟩ = [27].The corresponding creation and annihilation operators follow as where S is the overlap matrix.For brevity, the index refers to a multi-index, identifying the atomic position and orbital, e.g., = ( ⃗ R l , ⃗ R s , ) corresponds to lattice l, sublattice s and atomic orbital .Using the local basis set, we construct field operators, where the atomic orbitals,   (⃗ r − ⃗ R ls ) , are located at the equilibrium positions of the molecule or lattice.We build the usual Hamiltonian operator in second quantization, The transformation to the Bloch eigenstates, labeled by ⃗ k-vectors and band indices m, employs the usual transformation, where N c denotes the number of cells.The eigenstates U( ⃗ k) satisfy the secular equations, . Translational invariance of the Hamiltonian and the overlap matrices, e.g., H(l, l � ) = H(l � − l) , is used to Fourier transform and subsequently diagonalize the Hamiltonian, such that, To proceed with the derivation of the electron-phonon couplings we consider the perturbation of the Hamiltonian due to a deformation of the lattice.The corresponding Hamiltonian variation can be written as, where H(⃗ r) is the perturbation of the single-particle Hamil- tonian (2).Note that the field operators are defined in terms of the equilibrium positions and are thus not changing with atomic displacements.Therefore, the explicit form of the variation due to phonon displacements reads, where is the displacement operator in terms of the phonon states, labeled by the ⃗ q-vector and the polarization, .The mass of atom s is denoted by m s .The phase factors associated with the sublattice positions, ⃗ R s , are sometimes included in the polarization states, ẽ s (⃗ q) = e  s (⃗ q)e i⃗ q⋅ ⃗ R s , depending on the definition of the dynamical matrix [28].In phonopy [29], the phases are consistent with Eq. ( 8).
The last steps of our derivation involve writing the matrix elements in Eq. ( 7) in terms of variations of the usual TB matrix elements.This is accomplished by the identity, ( 6) which only shows sublattice indices in order to simplify the notation.In Eq. ( 9), we have explicitly used the two-center approximation of the DFTB matrix elements, allowing to remove the summation over s ′′ , and the relationship A consequence of this approximation is that the on-site ( s = s � ) electron-phonon coupling matrix- elements are zero.Further simplifications are possible by inserting the identity, Finally, it is convenient to rotate to the Bloch eigenstates (5) of the unperturbed Hamiltonian.Equation ( 10) contains four terms: the first and the third term lead to expressions containing, which, after the change of variables, �� = � − , can be transformed to In the second and the fourth term of Eq. ( 10), S −1 can be removed using the secular equation in favor of  n ( ⃗ k) and  m ( ⃗ k � ) .Collecting everything it is easy to obtain the equation of motion in the Heisenberg picture, where the electron-phonon couplings can be expressed as, (10) The difference of terms in ⃗ k and ⃗ k + ⃗ q within the curly braces is a consequence of the fact that the coupling depends on the relative displacement ⃗ u � s � − ⃗ u s .On the other hand, the matrices U( ⃗ k) and U( ⃗ k + ⃗ q) map local orbitals to eigenvec- tors, where the incoming electron with momentum ⃗ k is scat- tered into ⃗ k + ⃗ q .In fact, Umklapp processes come naturally into play in Eq. ( 14), since the lattice summation is always satisfied up to a reciprocal lattice vector,  ⃗ k � , ⃗ k+⃗ q+G .Consequently, momentum vectors are always to be understood as crystal quasi-momentum vectors.Equation ( 14) is in agreement with the result of [24], equation (2.31).We observe once more that the perturbation in Eqs. ( 6) and ( 7) is evaluated using the unperturbed basis set.Different expressions are obtained if the basis set is displaced together with the atomic oscillations.This, for instance, comes out automatically when considering some form of perturbations to Eq. ( 4) leading to variations of the Hamiltonian and overlap matrix elements, rather than matrix elements of the perturbed Hamiltonian, as pointed out in [24].Interestingly, if one starts from the Hamiltonian (4), expresses the equation of motion for bs and then considers perturbations to the matrix elements H + H and S + S , after some lengthy algebra arrives at a final equation for the couplings which is identical to Eq. ( 14) except that term.See Supporting Information for a more detailed discussion.

Implementation and results
In order to demonstrate our approach, we consider graphene as an example material for which we calculated the electron-phonon couplings, relaxation times, and transport properties.Our implementation, which uses phonopy [29] and dftb+ [11], is available under https:// github.com/ CoMeT 4MatS ci/ dftbe phy.For all DFTB calculations discussed in this section, we used the matsci-0-3 parametrization [30] without self-consistent charges.We also tested other parametrizations and the influence of self-consistent charges, but we found no significant differences (see Supporting Information).A 7 × 7 supercell was constructed from an optimized unit cell ( 10 −5 a.u.maximum force tolerance and 48 × 48 k-points).From the latter, we obtained the lattice constant a = 2.467 Å.

Electronic structure and phonons
For the supercell, the Hamiltonian and the overlap matrices, H , ( s, � s � ) and S , ( s, � s � ) , were obtained from dftb+ in real space.Those are Fourier transformed with respect to the cell indices to get H s,s � ( ⃗ k) and S s,s � ( ⃗ k) .Solving a generalized eigenvalue problem for each k-point yields the electronic band structure  n ( ⃗ k) and the respective eigenstates U s � m ( ⃗ k) as explained in Sect.1.2.The resulting band structure is shown in Fig. 1b.One can see the characteristic Dirac cone, ( ⃗ k) = ±�v F k , in the vicinity of the K-point.We find the Fermi velocity v F to be approximately 0.8 × 10 6 m/s in close agreement with previ- ously reported values [4,31].
To obtain the phonon dispersion and polarization vectors, i.e., frequencies   (⃗ q) and polarizations ẽ s (⃗ q) , we used pho- nopy with the same supercell.Due to the P6/mmm symmetry of graphene, only one displacement ( d = 0.005 Bohr mag- nitude) has to be considered.The dispersion is shown in Fig. 1c.Three of the six branches are acoustic modes with   (⃗ q) → 0 for ⃗ q → Γ .Since graphene is a 2D material, one of those modes is a bending mode (or ZA mode) with a quadratic dispersion [32].The other two acoustic modes are linear in k and correspond to the longitudinal (LA) and transversal (TA) acoustic modes.The respective sound velocities are found as c LA = 23.6 × 10 3 m/s and c TA = 13.6 × 10 3 m/s, which are close to the values reported before [4,31].

Electron-phonon coupling
In the next step, each atom in the supercell was displaced in each direction by = 0.001 Å and the real-space gradients of the Hamiltonian and overlap matrices were found from a finite-difference scheme.The gradients are also Fourier transformed yielding H s,s � ( ⃗ k)∕ ⃗ R s and S s,s � ( ⃗ k)∕ ⃗ R s .Together with the previously obtained output of phonopy the electron-phonon couplings g  nm ( ⃗ k, ⃗ q) were computed via Eq.( 14).
In the long wavelength limit and close to the K-points, the electron-phonon couplings for the TA and LA modes have a characteristic dependence on ⃗ q [4,31].In Fig. 2a, we show q as a function of ⃗ q for the conduction-band, pho- non modes = TA, LA and ⃗ k close to one of the K-points.T h e ch a r a c t e r i s t i c l e n g t h  ⃗ q i s g i ve n by where m C is the mass of a carbon atom.We obtain the expected symmetry, but compared to previously reported results in Refs.[4,31] our couplings seem to be smaller, although it should be noted that a quantitative comparison is not straight forward.One potential issue is that due to the different Fermi velocities the k-point ⃗ k for which the couplings are evaluated is different in different studies.The couplings for the ZA and ZO modes vanish.

Relaxation times
There are different approximations of the relaxation times which in general are calculated from the electron-phonon couplings.Specifically, we used the self-energy relaxationtime approximation (SERTA) of the scattering rate [1]: here Ω BZ is the area of the Brillouin zone.f 0 ( , , T) and n B ( , T) denote the Fermi function and the Bose-Einstein distribution, respectively, which characterize the occupation of electronic states and the phonons in equilibrium.They depend on temperature T and chemical potential .In our case, the partial decay rates are computed using a Gaussian smearing function with constant width ( = 3 meV) instead of the -functions.The integral over the Brillouin zone in Eq. ( 15) is replaced by a sum over q-points.We use a q-mesh ( 15) with 200 × 200 Monkhorst-Pack points which are sampled around the Γ-point up to q max ≈ 0.0665 × 2 ∕a .This implies that we consider only intravalley scattering.Intervalley scattering has to be treated separately using a shifted q-mesh [4].
Figure 3 shows the scattering rates (i.e., inverse lifetimes) of the LA and TA modes as a function of electronic energy in the conduction-band for T = 300 K and = 100 meV.As expected for the high-temperature regime and for chemical potentials close to the Dirac point, the scattering rates of the acoustic modes are linear in the electronic energy [31,33].The optical modes become relevant for energies larger than ℏ LO,TO (Γ) .When the energy level relative to the conduction-band (  ⃗ k − E C ) is 0.4 eV, the magnitudes of the scattering times are larger by about 50% compared to earlier calculations reported in [4].This observation does not take into account that the Fermi velocity obtained in the DFTB calculations is smaller by ≈ 12%.

Transport properties
Having the relaxation times, the conductivity and mobility of electrons and holes can be calculated by solving the Boltzmann transport equation [e.g.[1]].The conductivity tensor is given by ( 16) here A uc is the area of the unit cell and v n ( ⃗ k) and v n ( ⃗ k) are the charge carrier group velocities in the directions and ; for instance, v n ( ⃗ k) =  n ( ⃗ k)∕�k  .The excess chargecarrier densities n el and n h at a given temperature T and a chemical potential can be defined as ( c = el, h) Finally, the mobility is given by the ratio of conductivity and density, In our calculations, the integrals over ⃗ k are replaced by sums over ⃗ k-points.We used the same q-mesh for the calculation of the relaxation times described above and a k-mesh with 400 × 400 Monkhorst-Pack points.For the latter, only irreducible k-points were considered.The electron-phonon couplings were calculated for energies Consequently, the transport properties presented here do not include intervalley scattering.
Figure 4a shows the carrier concentration as a function of the chemical potential at T = 300 K. Since the consid- ered chemical potentials are still close to the Dirac point, the results can to a large extent be described using Dirac electrons.To illustrate this, the behavior resulting from the linear dispersion of the electronic bands around the K-point and of the acoustic phonon modes is included in Fig. 4 [31,33].This also applies to the conductivity, which is shown in Fig. 4b.Using the constant relaxation-time approximation (CRTA) the conductivity is seen to monotonously increase with the chemical potential, while it saturates within the SERTA.This stark difference is plausible when considering the linear dependence on energy of the scattering rate as discussed before.For the analytic SERTA behavior, we obtained −1 ( ) from Fig. 3. Combin- ing the carrier density and the conductivity, one obtains the mobility as shown in Fig. 4c.In the SERTA at 300 K, the mobility decreases approximately as 1/n [33].At n = 10 12 cm −2 we obtain for the mobility ≈ 1.3 × 10 5 cm 2 ∕ (Vs) which is similar to previously reported values [4].It should be noted that the inclusion of intervalley scattering leads to a mobility which is reduced by a factor of two compared to the pure intravalley case for a reasonable carrier density [4].For completeness, Fig. 4d shows the mobility as a function of temperature at fixed density n ≈ 4 × 10 12 cm −2 .( 17)

Conclusions
In summary, we presented a derivation of electron-phonon couplings within the general non-orthogonal tight-binding Hamiltonian.The final expression, Eq. ( 14), is in agreement with previously found results.However, depending on the choice of the basis set, i.e., centered at equilibrium positions or co-moving with the displaced atoms, one obtains slightly different expressions.The impact of this choice will be investigated elsewhere.The derivation is also discussed in the context of DFTB which allows for addressing large systems with reasonable accuracy.
We implemented the calculations of electron-phonon couplings in a Python package, dftbephy, which also provides the functionality to obtain relaxation times and transport properties within the Boltzmann transport framework.As an example, we show results for graphene.We find an overall good agreement with previously reported results.Using different parametrizations and self-consistent charges only has a minor influence on band-structures, phonon dispersions and relaxation times.For more complex materials the impact of SCC might be more pronounced, in particular if different elements are present.Overall, for materials with a large number of atoms in the unit cell, using DFTB to calculate electron-phonon couplings and transport-related properties provides an efficient approach to evaluate the performance of those materials for applications in electronic devices.

Fig. 1 a
Fig. 1 a Optimized geometry of graphene as a 7 × 7 supercell.Blue rhombus denotes the central unit cell.b The electronic band structure calculated with DFTB.Fermi level E F is shown with dashed line.c The phonon dispersion along Γ-M-K-Γ

Fig. 3
Fig. 3 Inverse lifetimes at 300 K as a function of energy k relative to the conduction-band minimum E C .Fermi level is 100 meV above E C

Fig. 4
Fig. 4 Numerical results (symbols) and analytic behavior (lines): a Carrier concentration as a function of the chemical potential for graphene.b Conductivity = xx + yy as a function of the chemi-