DFT2kp: effective kp models from ab-initio data

The $\mathbf{k}\cdot\mathbf{p}$ method, combined with group theory, is an efficient approach to obtain the low energy effective Hamiltonians of crystalline materials. Although the Hamiltonian coefficients are written as matrix elements of the generalized momentum operator $\mathbf{\pi}=\mathbf{p}+\mathbf{p}_{{\rm SOC}}$ (including spin-orbit coupling corrections), their numerical values must be determined from outside sources, such as experiments or ab initio methods. Here, we develop a code to explicitly calculate the Kane (linear in crystal momentum) and Luttinger (quadratic in crystal momentum) parameters of $\mathbf{k}\cdot\mathbf{p}$ effective Hamiltonians directly from ab initio wavefunctions provided by Quantum ESPRESSO. Additionally, the code analyzes the symmetry transformations of the wavefunctions to optimize the final Hamiltonian. This is an optional step in the code, where it numerically finds the unitary transformation $U$ that rotates the basis towards an optimal symmetry-adapted representation informed by the user. Throughout the paper, we present the methodology in detail and illustrate the capabilities of the code applying it to a selection of relevant materials. Particularly, we show a"hands-on"example of how to run the code for graphene (with and without spin-orbit coupling). The code is open source and available at https://gitlab.com/dft2kp/dft2kp.


I. INTRODUCTION
The band structure of crystalline materials defines most of its electronic properties, and its accurate description is essential to the development of novel devices.For this reason, the ab initio density functional theory (DFT) [1,2] provides one of the most successful tools for the development of electronics, spintronics, optoelectronics, etc.The DFT methods have been implemented in a series of codes (e.g., Quantum ESPRESSO [3,4], VASP [5], Wien2K [6], Gaussian [7], DFTB+ [8], Siesta [9,10], ...), which differ by the choice of basis functions (e.g., localized orbitals or plane-waves), pseudo-potential approximations, and other functionalities.Nevertheless, all DFT implementations provide methods to obtain the equilibrium (relaxed) crystalline structure, phonon dispersion, and electronic band structures.Complementary, few bands effective models are essential to further study transport, optical, and magnetic properties of crystalline materials.These can be developed either via the tightbinding (TB) [11][12][13] or k •p method [14,15], which complement each other.
On the one hand, the TB method has an "atomistic" nature, since it is built upon localized basis sets (e.g., maximally-localized Wannier functions [16], or atomic orbitals), which makes this method optimal for numerical modeling of transport, optical and other properties of complex nanomaterials [17][18][19][20].
On the other hand, the k • p method uses basis sets of extended waves, which are exact solutions of the Hamiltonian at a quasi-momentum of interest, typically at a high symmetry point of the Brillouin zone.While this characteristic may limit the k • p description to a narrow region of the energy-momentum space, the k • p Hamiltonians are easier to handle analytically and, especially, are very suitable to study mesoscopic systems using the envelope function approximation [21][22][23][24][25][26].For example, the k.p framework has been successfully applied to study nanostructures (quantum wells, wires, and dots) [27][28][29], topological insulators [30][31][32], spin-lasers [33,34], polytypism [35][36][37], as well as a large variety of two-dimensional van der Waals materials [38][39][40][41].Moreover, recent developments in the field of transition metal dichalcogenides (TMDCs) have combined DFT and k•p methodologies to explore the valley Zeeman physics in TMDC monolayers and their van der Waals heterostructures [42][43][44][45].
Both the TB and k • p Hamiltonians are defined in terms of arbitrary coefficients.In the TB case, these are local site energies and hopping amplitudes described by Slater-Koster matrix elements [11].For the k • p Hamiltonians, these are the Kane [46,47] and Luttinger [48] parameters, which are matrix elements of the momentum and spin-orbit coupling operators.In both methods (TB or k •p), the values of these arbitrary coefficients must be determined from outside sources, which strongly depend on the size and analytical properties of the particular model Hamiltonian.For instance, early studies within the k • p framework have shown that for parabolic single band descriptions, or weakly coupled models, it is possible to write the quadratic coefficients in terms of effective masses, which can be experimentally determined by cyclotron resonance experiments [46,[49][50][51][52].Moreover, energy splittings, such as band gaps, can be directly determined from optical experiments [53][54][55][56][57].For III-V semiconductors with zinc-blend structure and nitride-based wurtzite compounds, a useful database for k • p parameters inspired by experimentally available datasets can be found in Ref. [58].Conversely, for k • p Hamiltonians that do not allow analytical solutions, but still have a low number of bands (∼ 10), it is possible to perform numerical fitting techniques to DFT calculations [39,41,[59][60][61][62][63][64][65].For larger k • p Hamiltonians (> 30 bands), fitting proce-dures may also be applied [66,67] or directly extracted from first principles calculations, since the only matrix elements involved are linear in momentum [68][69][70].Interestingly, these large band k • p models can even be used to supplement and speed up first principles calculations, as demonstrated in Refs.[68][69][70].In TB models, fitting procedures can also be applied to obtain the unknown parameters [71][72][73][74][75]. Conversely, fully automated procedures, integrated within ab initio codes, such as the wannier90 code [76,77], use localized Wannier functions computed from the DFT wave functions to calculate TB parameters.Moreover, explicit calculations of the Slater-Koster matrix elements are implemented in the paoflow [78] and DFTB+ [8] codes.
While it is possible to extract k • p models from a Taylor expansion on top of a TB model (e.g., via the code tbmodels [79]), there are no versatile implementations to calculate the k • p Kane (linear in k) and Luttinger (quadratic in k) parameters directly from the DFT wavefunctions [80].To calculate the k•p matrix elements from the DFT wavefunctions, one needs to account for how the wavefunctions are represented in the DFT code [69,81].For instance, Quantum ESPRESSO and VASP implement pseudopotential approximations within the Projector Augmented Wave (PAW) method [82][83][84][85].Fortunately, Quantum ESPRESSO already provides a routine to calculate matrix elements of the velocity operator (which is sufficient to obtain k • p models, as we see in Section B).Indeed, recently, Jocić and collaborators [86] have successfully calculated k • p models directly from QE's wavefunctions (see disclaimer at our Conclusions).
In this paper, we present an open-source code that automatically calculates the numerical values for the k • p Kane and Luttinger parameters using the wavefunctions provided by Quantum ESPRESSO (QE).For this purpose, first, we develop a patch to instruct QE to calculate and store the matrix elements of the generalized momentum π = p + p SOC , which includes the spin-orbit corrections.Together with the eigenenergies E 0 n at k 0 , the matrix elements of π for a selected set of N bands define the effective k • p Hamiltonian H N ×N (k) for k near k 0 .Our python package reads these matrix elements and QE's wavefunctions |n⟩ to automatically build H N ×N (k) using Löwdin's partitioning [87] for the folding down of all QE bands into the selected N bands subspace.Additionally, the user has the option to improve the appearance (or form) of the effective Hamiltonian via a symmetry optimization process aided by the qsymm package [88], which builds the symbolic Hamiltonian via group theory and the method of invariants.To illustrate the capabilities of our code, we show here a step-by-step "hands-on" tutorial on how to run the code for graphene, and later we present results for selected materials [zincblende, wurtzite, rocksalt, transition metal dichalcogenides (TMDC), and others].In all cases, the modeled band structure matches remarkably well the DFT data at low energies near the expansion point k 0 .Our code is open source and available at the gitlab repository [89].This paper is organized as follows.In Section II we present our methodology starting with a brief review of the k • p method, Löwdin partitioning, the method of invariants, the symmetry optimization process, and the calculation of matrix elements using the DFT data.Next, in Section III, we show the code in detail using graphene as a practical example.Later, in Section IV, we illustrate the results of the code for zincblend (GaAs, CdTe, HgTe), wurtzite (GaP, GaN, InP), rock-salt (SnTe, PbSe), a TMDC (MoS 2 ), and other materials (Bi 2 Se 3 , GaBiCl 2 ).We finish the paper with an overview of the results in Section V, and the conclusions.

II. METHODS
Our goal is to obtain the numerical values for the coefficients of k • p effective Hamiltonians [14,15].Namely, these are the Kane [46,47] and Luttinger [48] parameters.To present our approach to this calculation, let us start by briefly describing its fundamental steps.First, we review the k • p method to show that these coefficients depend only upon matrix elements of the type P m,n = ⟨m| π |n⟩, where π = p + p SOC is the generalized momentum operator with the spin-orbit corrections, and {|n⟩} is the set of numerical wavefunctions obtained from the ab initio DFT simulations (e.g., via Quantum ESPRESSO [3,4]).However, the numerical DFT basis given by {|n⟩} does not match, a priori, the optimal symmetry-adapted basis set that yields the desired form for the effective k • p Hamiltonian.Therefore, to properly identify the Kane and Luttinger parameters, we perform a symmetry optimization, which rotates the arbitrary numerical basis into the optimal symmetry-adapted form.This symmetry optimization is performed via group theory [90,91] by enforcing that the numerical DFT basis transforms under the same representation of an optimal symmetry-adapted basis, which is informed by the user.
In summary, the algorithm steps are: 1. Read the QE/DFT data: energies E 0 n and eigenstates |n⟩ at the selected k 0 point.

Calculate or read the matrix elements of P m,n =
⟨m| π |n⟩ for all bands (m, n).

Select the bands of interest (set A)
. The code will identify the irreducible representations of the bands using the IrRep python package [92], and present it as a report to the user.Additionally, the code calculates the model folded down into the selected set A via Löwdin partitioning.

4.
Build the optimal effective model from symmetry constraints using the Qsymm python package [88] under an optimal symmetry-adapted basis informed by the user.This optimal basis must be in a set of representations equivalent to the ones identified in Step 3.

5.
Calculate the representation matrices for the symmetry operators in the original QE basis |n⟩.The code verifies if the representations of the numerical QE basis are equivalent to the representations of the optimal symmetry-adapted basis from step 4.

6.
Calculates the transformation matrix U that rotates the original QE basis into the optimal symmetry-adapted basis set in step 4. Applies the transformation U and calculates the optimal symmetry-adapted numerical effective Hamiltonian.
7. Convert values from Rydberg atomic units into meV and nm units, and present a report with values for the k • p parameters.
In the next sections, we describe the relevant details of the steps above, but not following the algorithmic order above.More specifically, in Section II A, we briefly review the k • p formalism to show that P m,n = ⟨m| π |n⟩ plays a central role in our approach.Incidentally, we introduce the folding down via Löwdin partitioning [87].Next, we define what is the optimal symmetry-adapted form of the Hamiltonian via the method of invariants [15,93] in Section II B. In Section II C, we present the symmetry optimization approach to calculate the transformation matrix U that yields our final Throughout the paper we use atomic Rydberg units (a.u.), thus the reduced Planck constant, bare electron mass and charge are ℏ = 2m 0 = e 2 /2 = 1, the permittivity of vacuum is 4πε 0 = 1, the speed of light is c = 2/α ≈ 274, and α ≈ 1/137 is the fine structure constant.
We are interested in the effective Hamiltonian near a high-symmetry point k 0 of the Brillouin zone.Therefore, we write the quasi-momentum as κ = k 0 + k, such that k is the deviation from k 0 .The Bloch theorem allow us to decompose the wavefunction as ψ κ (r) = e ik•r ϕ k0,k (r), with ϕ k0,k (r) = e ik0•r u k0+k (r), where u k0+k (r) ≡ u κ (r) is the periodic part of the Bloch function, while ϕ k0,k (r) carries the phase given by k 0 and obeys the Schrödinger equation where H 0 is the Hamiltonian at k = 0, V (r) is the periodic potential, H ′ (k) carries the k-dependent contributions that will be considered as a perturbation hereafter, π is the generalized momentum that includes the spin-orbit contributions (SOC), and σ = (σ x , σ y , σ z ) are the Pauli matrices for the electron spin.For simplicity, we consider only leading order corrections of the fine structure terms.Namely, at k = 0, the H SR carries the scalar relativistic terms, composed by the Darwin, H D = α 2 8 ∇ 2 V (r), and the mass-velocity corrections, H MV = −α 2 p 4 /4.In the ab initio DFT data, these are implied in the numerical eigenvalues E 0 n of H 0 .For finite k ̸ = 0, we keep only the SOC contribution in π, and neglect the higher order mass-velocity corrections (see Appendix A).
The DFT data, as shown in the next section, provide us with a set {|n⟩} of eigenstates of H 0 , i.e.H 0 |n⟩ = E 0 n |n⟩.From this crude DFT basis, we define an all bands model H DFT all (k), with matrix elements where P m,n = ⟨m| π |n⟩.We refer to this as the crude model because it is calculated from the original numerical DFT wavefunctions, which do not have an optimal symmetry-adapted form (more detail in Section II C).Nevertheless, it already shows that E 0 n and P m,n are central quantities, and both can be extracted from DFT simulations, as shown in Section II D.
Next, we want to fold down H DFT all into a subspace of N bands near the Fermi energy to obtain our reduced, but still crude, effective model H DFT N ×N .This is done via Löwdin partitioning [15,87,93].First, the user must inform the set of N bands of interest, which we refer to as set A. Complementary, the remaining remote bands compose the set B. Considering the diagonal basis H 0 |n⟩ = E 0 n |n⟩, and the perturbation H ′ (k), the Löwdin partitioning leads to the effective Hamiltonian H DFT N ×N defined by the expansion with Here, the indices m, n ∈ A run over the bands we want to model (set A), while r ∈ B run over the remote bands.The expansion above is shown up to second order in H ′ , but higher order terms can be found in Ref. [15].Alternatively, the recent python package pymablock [94] implements an efficient numerical method to compute the Löwdin partitioning to arbitrary order.
B. The optimal symmetry-adapted form of H The selection rules from group theory allow us to identify which matrix elements of an effective Hamiltonian are finite [90].More interestingly, the method of invariants [15,93] can be used to directly obtain the most general form of H optimal N ×N (k) allowed by symmetry.To define this form, consider a Taylor series expansion where h i,j,l are constant matrices that multiply the powers of k = (k x , k y , k z ) as indicated by its indices i, j, l = {0, 1, 2, . . .}.To find the symmetry allowed h i,j,l , we recall that the space group G of the crystal is defined by symmetry operations that keep the crystalline structure invariant.Particularly, at a high symmetry point κ = k 0 , one must consider the little group G k0 ∈ G of symmetry operations that maintain k 0 invariant (the star of k 0 ).Hence, H optimal N ×N (k) must commute with the symmetry operations of G k0 .Namely, where D ψ (S) are the representation matrices for each symmetry operator S ∈ G k0 in the subspace defined by the wavefunctions of set A, and D k (S) are the representation matrices acting on the vector k = (k x , k y , k z ).The set of equations defined by this relation for all S ∈ G k0 leads to a linear system of equations that constrain the symmetry allowed form of H optimal N ×N (k), i.e., it defines which of constant matrices h i,j,l are allowed up to a multiplicative factor.Ultimately, these multiplicative factors are the Kane and Luttinger parameters that we want to calculate numerically.
The python package Qsymm [88] implements an efficient algorithm to find the form of H optimal N ×N (k) solving the equation above and returns the symmetry allowed h i,j,l .Qsymm refers to these as the Hamiltonian family.To perform the calculation, the user must inform the representation matrices D ψ (S) for the generators of G k0 .Notice that the choice of representation is arbitrary, and different choices lead to effective Hamiltonians with different forms.This ambiguity is the reason the next step, symmetry optimization, is necessary.

C. Symmetry optimization
In the previous section, the matrix representations for generators S ∈ G k0 are implicitly written in an optimal symmetry-adapted basis, which we will now label with an O index, as in {|n O ⟩}, to distinguish from the crude DFT numerical basis, which we now label with an C index, as in {|n C ⟩}.The matrix representations of S written in these two bases are equivalent up to a unitary transformation U , i.e.D O (S) = U •D C (S)•U † .Indeed, this same matrix U transforms the crude DFT numerical Hamiltonian into the desired optimal symmetry-adapted form, i.e.
Therefore, our goal here is to find this transformation matrix U .
For each symmetry operator S i ∈ G k0 , let us define C i ≡ D C (S i ) and O i ≡ D O (S i ) as the representation matrices under the original numerical DFT basis (C), and under the desired optimal symmetry-adapted representation (O), respectively.For irreducible representations, this U is unique (modulo a phase factor) and an efficient method to obtain it was recently developed [95] and used in Ref. [86] to transform the effective model into the desired form.The procedure described in Ref. [95] is exact but relies on a critical step where one has to find for which indices (a, b) the weight matrix r a,b is finite.For transformations between irreps, any of the finite r a,b lead to equivalent unitary transformations.However, for transformations between reducible representations, one needs to identify, within the set of finite r a,b , the ones that yield nonequivalent transformation matrices that combine to form the final transformation matrices U .This can be a complicated numerical task.Here, instead, we propose an alternative method that applies more easily to reducible representations and allows us to obtain the transformation matrix U with a systematic approach.Next, we describe the method, and later in Sec.III C we illustrate its capabilities using the spinful graphene example.
The set of unitary transformations k0 compose a system of equations for U .These can be written in terms of its matrix elements in a linearized form that reads as Defining a vector where N is the order of the representations (number of bands in set A), allow us to cast the equation above as , and 1 N as the N × N identity matrix.Since the same similarity transformation U must apply for all S i , we stack each Q i into a rectangular matrix The full set of equations now read as is a linear combination of the nullspace {v j } of Q, with coefficients c j and nullity N Q .The matrix U can be recovered from the elements of V , which follow from its definition above.If u j is the matrix reconstructed form of v j , we can write U = N Q j=1 c j u j .Additionally, it is interesting to consider anti-unitary symmetries.These can be either the time-reversal symmetry (TRS) itself, or combinations of TRS and space group operations (magnetic symmetries) [90,91].For instance, in spinful graphene neither TRS nor spatial inversion are symmetries of the K point, but their composition is an important symmetry that enforces a constraint on the allowed SOC terms (see Sec. III C).Following a notation similar to the one above, let us refer to these magnetic symmetries as Ci = D C ( Si )K ≡ Ci K and Ōi = D O ( Si )K ≡ Õi K, where K is the complex conjugation, and ( Ci , Õi ) are the unitary parts of ( Ci , Ōi ).Now the basis transformation for these symmetries read as Õi = U * • Ci •U † , where we choose to apply K to the left (this choice is for compatibility with the python package IrRep [92]).To add this equation to the Q matrix above, we consider U and U * as independent variables.Then, as above, it follows the linearized form In all cases, the expression for the transformation matrix is U = N Q j=1 c j u j , where the coefficients c j are so far undefined.To find these coefficients c j , we numerically minimize the residues R({c j }) The global minima of these residues, R({c j }) = R({c j }) ≡ 0, yields a solution U ({c j }), such that small perturbations to the coefficients c j → c j + δc j lead to quadratic deviations from the minima, e.g., R ∝ |δc j | 2 .This procedure opens a question of whether the solution U ({c j }) at the global minima is unique.
Since U represents a transformation between two basis sets (e.g., |n O ⟩ = U |n C ⟩), it expected to be unique.However, the problem here is formulated such that we explicitly have the eigenstates |n C ⟩ that compose the crude DFT basis set C, while for the optimal symmetry-adapted basis set O we know only how we expect the eigenstates |n O ⟩ to transform under the symmetry operations of the group.Therefore, instead of solving for U directly from the linear basis transformation |n O ⟩ = U |n C ⟩, we rely on the quadratic equations for the transformation between the symmetry operators (e.g., D O (S) = U •D C (S)•U † ), or their linearized forms in Eq. ( 8) and Eq. ( 9).First, consider that O and C refer to distinct, but equivalent irreps.As emphasized in [95], it follows from Schur's lemma that the transformation U is unique modulo a phase.Indeed, for the unitary constraints, Next, without loss of generality, let us consider that O and C refer to reducible representations already cast in block-diagonal forms.In this case, the solution U = U 1 ⊕ U 2 ⊕ • • • also takes a block-diagonal form, where each block U j corresponds to a transformation within a single irrep subspace.It follows that each U j is unique modulo the phases above.The overall global phase of U does not affect the calculation of our matrix elements.However, the arbitrary relative phases between the blocks U j might lead to illdefined phases of matrix elements between eigenstates of different irreps if the anti-unitary symmetries are not informed.In contrast, if anti-unitary symmetries are used, the undefined phase factor in the matrix elements is just a sign.

D. Matrix elements via DFT
As shown above, our approach to obtain a k • p model directly from the DFT data relies on two quantities: (i) the band energies E 0 n at the k • p expansion point k 0 ; and (ii) the matrix elements P m,n = ⟨m| π |n⟩ also calculated at k 0 for all bands {|n⟩}.The band energies E 0 n are a straightforward output of any DFT code.Therefore, here we discuss only the calculation of P m,n = ⟨m| π |n⟩.
We focus on the Quantum ESPRESSO (QE) [3,4] implementation of ab initio DFT [1,2].There, the Hamiltonian is split into the core and intercore regions via the Projector Augmented Wave (PAW) method [82][83][84], which is backward compatible with ultrasoft (USPPs) [83,96] and norm-conserving pseudo-potentials (NCPP) [97][98][99].In these approaches, the atomic core region is replaced by pseudopotentials, which are constructed from single-atom DFT simulations with the Dirac equation in the scalar relativistic or full relativistic approaches.Thus, for molecules or crystals, QE solves a pseudo-Schrödinger equation, with the atomic potentials replaced by the pseudopotentials.Here we shall not go through the details of the PAW and pseudopotential methods.For the interested reader, we suggest Refs.[82][83][84].Instead, for now, it is sufficient to conceptually understand that QE provides numerical solutions for the Schrödinger equation with the fine structure corrections, which can be expressed by the Hamiltonian where H SR = H D + H MV contain the Darwin and massvelocity contributions, as presented above, and the last term is the spin-orbit coupling.

Matrix elements of the velocity
Fortunately, the QE code already provides tools to calculate the matrix elements of the velocity operator where we neglect the mass velocity corrections (see Appendix A).Thus, we find that P m,n = ⟨m| π |n⟩ ≈ ⟨m| 1 2 v |n⟩.The calculation of P m,n is already partially included in the post-processing tool bands.x(file PP/src/bands.f90),within the write_p_avg subroutine (file PP/src/write_p_avg.f90).This calculation includes the necessary PAW, USPPs, or NCPPs corrections, which are critical for materials where the wavefunction strongly oscillates near the atomic cores [100].However, the write_p_avg subroutine only calculates |P m,n | 2 for m in the valence bands (below the Fermi level) and n in the conduction bands (above the Fermi level).To overcome this limitation, we have built a patch that modifies bands.f90and write_p_avg.f90to calculate P m,n for all bands.This leads to a modified bands.xwith options to follow with its original behavior or to calculate P m,n according to our needs.This is controlled by a new flag lpall = False/True added to the input file of bands.x in addition to the lp = True.Its default value (lpall = False) runs bands.xwith its original code, while the option lpall = True instructs bands.x to store all P m,n into the file indicated by the input parameter filp.
In general, it is preferable to patch QE to use the full P m,n , since the calculation is faster and more precise.Nevertheless, if the user prefers not to apply our patch to modify QE, our code can calculate an approximate P m,n using only the plane-wave components outputted by the QE code.In this case, we consider that the pseudo-wavefunction is a reasonable approximation for the all-electron wavefunction, thus neglecting PAW corrections, which are necessary to account for SOC.Therefore, under this approximation, P m,n ≈ ⟨m| p |n⟩.The relevance of these PAW/SOC corrections to P m,n are presented in the example shown in Sec.IV B 1. Within this approximation, the wavefunction ψ n,k (r) for the band n at quasi-momentum k, and P m,n read as where c n (G) are the plane-wave expansion coefficients (spinors in the spinful case), Ω is the normalization volume, and G are the lattice vectors in reciprocal space.To implement this calculation, and the one shown next, we use the IrRep python package [92], since it already has efficient routines to read and manipulate the QE data.

Matrix elements of the symmetry operators
To calculate the matrix elements of the symmetry operators, it is sufficient to consider ψ n,k (r) from Eq. (12).In this case, it is safe to neglect PAW corrections, since they must transform identically to the plane-wave parts under the symmetry operations of the crystal space group.For a generic symmetry operation S ∈ G k0 , its matrix elements read as Using the plane-wave orthogonality, one gets where S −1 is the inverse of S, and

III. HANDS-ON EXAMPLE: GRAPHENE
In this section, we present a detailed example and results for spinless graphene, and a shorter discussion on spinful graphene in Sec.III C to illustrate the case of transformations between reducible representations.Graphene [101,102] is nowadays one of the most studied materials due to the discovery of its Dirac-like effective low energy model, which reads as H = ℏv F σ•k. Here, the σ Pauli matrices act on the orbital pseudo-spin subspace, k = (k x , k y ) is the quasi-momentum, and v F is the Fermi velocity, which is the unknown coefficient that we want to calculate in this example.For this purpose, we follow a pedagogical route in this first example.First, we present the symmetry characteristics of the graphene lattice and its wavefunctions at the K point.Then, we show the results for the representation matrices and Hamiltonian in the crude and optimal symmetry-adapted basis to illustrate how the symmetry optimization of Section II C is used to build the optimal symmetry-adapted Hamiltonians and identify the numerical values for its coefficients.Later, in Section III B we show a step-by-step tutorial on how to run the code.This example was chosen for its simplicity, which allows for a clear discussion of each step.Later, in Section IV we present a summary of examples for other materials of current interest.
Before discussing the details, we summarize the results for the band structure of graphene in Fig. 1, which compares the DFT data with our two main models.The black lines are calculated from the all bands model from Eq. ( 4), which uses the matrix elements P m,n in the original crude DFT basis without further processing.In contrast, the red lines are the band structure calculated with the folded-down Hamiltonian for a set A composed by the two bands near the Fermi energy that defines the Dirac cone, and considers the symmetry optimization process to properly identify the k • p parameters.This optimal symmetry-adapted Hamiltonian is shown in Eq. ( 21) below, and the numerical value for its parameters is shown at Step 7 in Section III B.

A. Overview of the theory and symmetry optimization
The crystal structure of graphene is a hexagonal monolayer of carbon atoms, as shown in Figs.1(a) and 1(b), which is invariant under the P6/mmm space group (#191).However, since its Dirac cone is composed of p z orbitals only, it is sufficient to consider the C 6V factor group to describe the lattice.Particularly, at the K point [see Fig. 1(c)], the star of K corresponds to the little group C 3V , which is generated by a 3-fold rotation C 3 (z) and a mirror M y .The Dirac bands of graphene are  4)] (black lines), and optimal symmetryadapted model [Eq.( 21)] for the two bands forming the Dirac cone (red).Here, the QE/DFT simulation was performed with 300 bands.
To build the optimal symmetry-adapted effective model via the method of invariants, we need to specify a basis and calculate the matrix representation of the symmetry operations mentioned above.Since the wavefunctions of the Dirac cone transform as the irrep E of C 3V , a naive choice would be A unconv ={|XZ⟩ ,|Y Z⟩}, which corresponds to a set A in Section II A. This choice of basis refers to a possible C representation in Section II C, and it yields where θ = 2π/3.Here H unconv is obtained via Qsymm up to linear order in k, for brevity.While the eigenenergies of H unconv represent correctly the Dirac cone as , the Hamiltonian H unconv takes an undesirable unconventional form.
A more convenient choice is A conv = {|(X + iY )Z⟩ ,|(X − iY )Z⟩}, which is illustrated in Figs.1(a) and 1(b).This choice of basis leads to where k ± = k x ±ik y .Now, up to linear order in k, we see that H conv ≈ c 0 +c 1 σ •k, where σ act on the subspace set by A conv , and we identify c 1 = ℏv F .Additionally, the kquadratic terms that lead to trigonal warping corrections.Notice that both choices, A unconv and A conv , are equivalent representations, but the conventional one leads to the familiar form of the graphene Hamiltonian.These two basis sets are related by an unitary transformation U , such that Next, let us analyze the set A QE of numerical wavefunctions from QE. Do they correspond to A QE = A conv or A QE = A unconv ?The answer is neither.Since it is a raw numerical calculation, typically diagonalized via the Davidson algorithm [104], a degenerate or nearly degenerate set of eigenstates might be in any linear combination of its representative basis.Therefore, the symmetry optimization step is essential to find the matrix transformation U that yields A conv = U • A QE .To visualize this, let us check the matrix representations of the symmetry operators above, and the effective Hamiltonian calculated from the crude QE data.For the symmetry operators, we find While this cumbersome numerical representation does not resemble neither A conv nor A unconv , our symmetry optimization process correctly finds a transformation matrix U that returns Finally, for the Hamiltonian, up to linear order in k and in the original QE basis, we find which takes a cumbersome form in this raw numerical basis.However, applying the transformation U , the symmetry adapted model becomes Here we identify ℏv F = 0.72 in Rydberg units, yielding v F = 0.83 × 10 6 m/s.The resulting band structure calculated from H optimal , including the k-quadratic terms, is shown as red lines in Fig. 1(d) and it matches well the QE/DFT data near K.

B. Running the code
The example presented here is available in the Examples/graphene-nosoc.ipynb notebook in the code repository, and shown in Algorithm 1.
Here we show only the minimal procedure to read the DFT data, build an effective model from the symmetry constraints, and calculate the numerical values for the model parameters.Complementary, the full code in Examples/graphene-nosoc.ipynb shows how to plot the data presented in our figures.
For now, we assume that the DFT simulation was successful.The suggested steps to run QE and prepare the data for our code is to run the calculation='scf' and calculation='bands' with pw.x.
Then, run bands.xto extract the bands from QE's output and store it in gnuplot format to plot the figures.Here, for graphene, we assume that the bands calculation was run for a path Γ − K − M with 30 points between each section, such that K is the 31st point in the list.
Next, we describe each step shown in Algorithm 1.
Step 1.After running QE, the first step is to read the DFT data from the QE's output folder.The command dft2kp.irrep(...) uses the python package IrRep [92] to read the data for the selected k point to be used in the k • p expansion, as indicated by the parameters kpt and kname.The data is read from the folder indicated by the parameter dftdir, while outdir and prefix refer to values used in the input file of QE's pw.x calculation.Additionally, the command dft2kp.irrep(...) also accepts extra parameters from the package IrRep (see code documentation).
Step 2. In step 2, the code will either read or calculate the matrix elements P m,n to build the effective models.If the user runs QE modified by our patch, the QE tool bands.xwill generate a file kp.dat that already contains the values for P m,n .In this case, the user must inform the name of this file via the parameter qekp.Otherwise, if qekp is omitted, our code calculates an approximate value for P m,n ≈ ⟨m| p |n⟩ from the pseudo-wavefunction of QE, as in Eq. ( 13), which neglects all SOC corrections.
Step 3. Next, the user must choose which set of bands will be considered to build the model.This is the set A in Section II A. In this example, we select bands 3 and 4, which correspond to the Dirac cone of graphene.The code analyzes the list of bands and identifies their irreducible representations (irreps) using the IrRep package [92].Here, the set A must contain only complete sets of irreps, otherwise the Löwdin perturbation theory would fail with divergences [see Eq. ( 5)], since the remote bands of set B would have at least one band degenerated with a band from set A. If this condition fails, the code stops with an error message.Otherwise, if set A is valid, the code outputs a report indicating the space group of the crystal (e.g., P6/mmm), the selected set of bands (e.g., [3,4]), their irrep (e.g., K 6 [103]), and degeneracy (2).

The report reads as
Space group 191: P6/mmm Verifying set A: [3 4] Band indices: [3,4] Irreps: (K6) Degeneracy: 2 Additionally, in this step, the code also calculates the crude effective model for the bands in set A via Löwdin partitioning [87].It stores the folded Hamiltonian in a Python dictionary (kp.Hdict) representing the matrices h i,j,l in the crude DFT basis that define H DFT (k) = i,j,l h i,j,k k i x k j y k l z .
Step 4. In step 4 we build the optimal symmetryadapted model using Qsymm [88], which solves Eq. ( 7) for the method of invariants.In Algorithm 1, we build the representations for the symmetry operations C 3 (z), M y , M z , and T I. Above we have discussed only the first two for simplicity.Here we also include the mirror M z , and the anti-unitary symmetry T I, which is composed of the product of time-reversal and spatial inversion symmetries.The mirror M z has a trivial representation D ψ (M z ) = −1, since the orbitals that compose the Dirac bands in graphene are all of Z-like (odd in z).The T I representation follows from A conv presented above by recalling that spinles time-reversal is simply the complex conjugation and the spatial inversion takes (X, Y, Z) → (−X, −Y, −Z).In this particular example, the T I symmetry does not play an important role, but it is essential for a spinful graphene example, as it constrains the SOC terms at finite k (see Sec. III C).The command dft2kp.qsymm(...) calls Qsymm to build the effective model from the list of symmetries, indicated by symm, up to order k 2 , as indicated by total_power.We recommend always using dim=3 [three dimensions for k = (k x , k y , k z )] because QE always work with the 3D space groups.Additionally, the command dft2kp.qsymm(...) accepts other parameters that are given to the Qsymm package (see code documentation).By default, this command outputs the optimal symmetry-adapted Hamiltonian, which matches the one in Eq. ( 21).
Step 5. Next, we start the symmetry optimization process.The first call kp.get_symm_matrices() calculates, via Eq.( 15), the matrix representation for all symmetry operators identified in the QE data by the IrRep package.However, neither QE nor IrRep account for the anti-unitary symmetries.Therefore, we call here the optional routine kp.add_antiunitary_symm(...), which manually adds the anti-unitary symmetry to the list of QE symmetries and matches it with the corresponding symmetry of Qsymm informed on its first parameter.In this example, we add the T I symmetry built with Qsymm above.This operator needs to be complemented with a possible non-symmorphic translation vector, which is zero in this case, as shown by the second parameter of kp.add_antiunitary_symm(...).Both calls, kp.get_symm_matrices() and kp.add_antiunitary_symm(...), calculate the matrix representations in the crude QE basis.
Step 6.To calculate the transformation matrix U , we compare the ideal matrix representations informed via Qsymm (object qs) and the crude QE matrix representations (object kp).
The call dft2kp.basis_transform(...) performs this comparison and returns an error if the symmetries in both objects do not match.More importantly, it calculates the transformation matrix U solving Eq. ( 8) and Eq. ( 9).The matrix U is stored in the object optimal.U.If the calculation of U is successful, the code applies U to rotate the h i,j,l terms in kp.Hdict from the crude DFT basis into the optimal symmetry-adapted basis.This allows for direct identification of the coefficients c n from Eq. ( 21), which are stored in optimal.coeffs.Additionally, the code builds the numerical optimal symmetry-adapted model and provides a callable object optimal.Heff(kx, ky, kz) that returns the numerical Hamiltonian H optimal N ×N for a given value of k = (k x , k y , k z ).
Step 7. At last, the code prints a report with the numerical values for the coefficients c n , which are summarized in Table I.As mentioned above, here we identify ℏv F = 0.72 a.u., yielding v F = 0.83 × 10 6 m/s after converting the units.
Algorithm 1 Minimal example for spinless graphene.C. Spinful graphene To complement the example above, we consider now the spinful graphene (full code available at Examples/graphene.ipynb [89]).In this case, due to the small spin-orbit coupling of graphene, the numerical DFT basis functions from QE mix two nearly degenerate irreps into an unintended reducible representation.Nevertheless, our symmetry optimization procedure can properly block diagonalize the symmetry operators according to the intended representation.To see this, let us first establish the ideal basis in proper ordering that leads to the block-diagonal form of the symmetry operators C 3 (z), M y , M z , and T I (considering the group generators only).Thus, considering the spin, the basis functions now read as {|(X + iY )Z, ↑⟩, |(X − iY )Z, ↓⟩, |(X − iY )Z, ↑⟩, |(X + iY )Z, ↓⟩}.Under the P6/mmm double space group [103,105], this set of basis functions transform as the sum of two bidimensional irreps [106], namely K7 ⊕ K9 .Under this basis, the symmetry operators listed above take a block-diagonal form, which are illustrated in the top row of Fig. 2. Algebraically, these read In contrast to the block diagonal form of the D ideal (• • • ) matrices above, the representation matrix for the C 3 (z) calculated with the crude DFT basis from QE takes the form Similarly, the crude DFT representation for M y , M z and T I also show non-block-diagonal forms in the central line of Fig. 2.
The algorithm described in Sec.II C builds a system of equations to find the transformation matrix U that yields D ideal (S) = U D QE (S)U † for all symmetry S of the group (i.e., S = {C 3 (z), M y , M z , T I} in this example).The Python code to implement this procedure is nearly identical to Algorithm 1, requiring only (i) the expansion of setA, in Step 3, to account for the 4 bands that compose the spinful Dirac cone (i.e., setA = [6, 7, 8, 9] in this Example); and (ii) the replacement of the symmetry matrices from Step 4 for the ones listed above.From these, in Step 6 we find the transformation matrix which precisely yields the transformation U D QE (S)U † = D optimal (S) ≡ D ideal (S), as illustrated in the bottom row of Fig. 2.
The model resulting from the considerations above read as where k 2 = k 2 x + k 2 y , k ± = k x ± ik y , and we omit k zdependent for 2D materials.Notice that if we do not consider the composed magnetic anti-unitary symmetry T I, the c 2 and c 5 terms above split into real and imaginary parts.Particularly for c 2 , the real part refers to matrix elements of p, while the imaginary part would carry contributions from p soc .Nevertheless, considering T I, these coefficients are expected to be real and the p soc contributions to the imaginary part vanish by symmetry.
The numerical values found for the parameters of H sfg in Eq. ( 34) are shown in Table II.The Fermi velocity matches the one from spinless graphene above, and we find that the intrinsic spin-orbit coupling is λ I = c 1 −c 0 ≈ 1 µeV, which is much smaller than its established value of λ I ≈ 24 µeV obtained via all-electron full-potential DFT implementations [107,108].This discrepancy is due to limitations of the pseudo-potentials used here with QE [109], which do not include d orbitals.Nevertheless, this example serves to show that, whenever two irreps are nearly degenerate, the DFT wavefunctions might always be mixed into reducible representations and the symmetry optimization procedure implemented here efficiently rotates the DFT basis back into ideal form that yields block-diagonal reducible representations.

IV. EXAMPLES
In this section, we briefly show the results for a series of selected materials without presenting a step-by-step tutorial as above.More details for each case below can be seen in the code repository.Here we consider examples of zincblende crystals (GaAs, HgTe, CdTe), wurtzite crystals (GaN, GaP, InP), rock-salt crystals (SnTe, PbSe), a transition metal dichalcogenide monolayer (MoS 2 ), 3D and 2D topological insulators (Bi 2 Se 3 , GaBiCl 2 ).Additional examples can be found in the code repository.In all cases, the resulting models agree well with the DFT bands near the k • p expansion point and low energies, as expected.The DFT parameters used in the simulations are presented in Appendix B.

A. Zincblende crystals
We consider well-known zincblende crystals: GaAs, CdTe and HgTe.These crystals are characterized by lattices that transform as the space group F 43m, but their low energy bandstructure concentrates near the Γ point, which can be described by the point group T d after factorizing the invariant subgroup of Bloch translations.The basis functions and effective Kane model for these materials are well described in the literature [14,15,91].
Here, let us simply summarize this characterization to establish a notation.
In all cases considered in this section, the first conduction band and the top valence bands transform either as S or P = (X, Y, Z) orbitals, and in terms of the crystallographic coordinates we define x ∥ [100], y ∥ [010], and z ∥ [001].In the single group T d , neglecting spin, the S-like orbitals transform accordingly to the trivial A 1 irrep of T d , while the P-like orbitals transform as the T 2 irrep.Including spin, the double group representation for the S-like orbitals become A 1 ⊗ D 1/2 = Γ6 , where D 1/2 is the spinor representation, and it yields the spin 1/2 basis functions |S ↑⟩ and |S ↓⟩.For the P-like bands one gets T 2 ⊗ D 1/2 = Γ8 ⊕ Γ7 , where Γ8 represents the basis functions of total angular momentum 3/2, and Γ7 has total angular momentum 1/2.These basis functions are listed in Table III.For GaAs and CdTe the conduction band is represented by Γ6 (S-type, and spin 1/2), the first valence band is composed of P-type orbitals with total angular momentum 3/2, which are described by the Γ8 irrep, and the split-off band contains P-type orbitals with total angular momentum 1/2, which defines the irrep Γ7 .In contrast, for HgTe the Γ6 and Γ8 are inverted due to fine structure corrections.
The basis from Table III diagonalizes the spinful effective Hamiltonian at k = 0, and leads to the well known extended Kane Hamiltonian [15].The expression for the 8 × 8 Hamiltonian H ZB is shown in Appendix C in terms of the coefficients c j following the output of the qsymm code, so that it matches Examples in our repository.There, the notation for the powers of k follows from Ref. [15], such that it can be directly compared to the extended Kane model shown in their Appendix C. The values for the coefficients c j are also shown in Appendix C. The band structures calculated from H ZB are shown in Fig. 3, which also shows the crystal lattice and the first Brillouin zone in Figs.3(a-b).In all cases, Figs.3(c-e), the blue dots represent the DFT results.The black lines are the crude model from Eq. 4, which includes all DFT bands and approaches a full zone description, but with a cost of a large N × N model with typical N ≫ 100.More importantly, the red lines represent effective 8 × 8 Kane model from H ZB , which matches well the DFT data at low energies and near Γ, as shown in the zoomed insets below each panel for GaAs [Fig.3(c)], HgTe [Fig.3(c)], and CdTe [Fig.3(c)].Particularly, for HgTe it is clear the band inversion between the Γ6 and Γ8 irreps.
Table III.Basis functions for zincblende crystals.The first column indicates the double group irreps for the T d point group at Γ, which are induced from the single group irreps in parenthesis.The second column lists the basis functions on the basis of total angular momentum, and the third column shows their expressions in terms of the symmetry orbitals (S, X, Y, Z) and spin (↑, ↓), which follows the definitions from Ref. [15].
The wurtzite crystals form a lattice that is characterized by the space group P6 3 mc, and the low energy band structure appears near the Γ point only.Near Γ, one can factorize the translations and the resulting factor group is the C 6V point group, which is generated by the C 6 rotation around the z-axis, and the mirror M x .Here, in terms of the crystallographic coordinates, x ∥ To illustrate the results for wurtzite materials, we consider the cases of GaN, GaP, and InP.Their band structures are shown in Figs.4(c-e).In all cases, the top valence bands are characterized by the irreps Here, A 1 is the trivial irrep of C 6V (single group), which represents S-like and Z-like orbitals, and E 1 is the vector representation of C 6V that contains (X, Y)-like orbitals.These are composed with the pure spinor representation D 1/2 to define the C 6V double group irreps Γ7 and Γ9 .Additionally, we consider two conduction bands, which are characterized by the irreps (A 1 + B 1 ) ⊗ D 1/2 = Γ8 ⊕ Γ9 .The orbital basis function for the B 1 irrep is odd under both C 6 and M x , its representation on group character tables is cumbersome, so one defines it as X(X 2 − 3Y 2 ) ≡ |V ⟩ [14].Ultimately, we consider the double group representations ordered as shown in Table IV.
Table IV.Basis functions for wurtzite crystals.The first column shows the double group irreps of C6V , which are induced from the single group irrep between parenthesis.The second column shows the basis representation in terms of the spherical harmonics Y m l and spin (↑, ↓), while the third column shows the representation in terms of the orbitals (S, X, Y, Z, V), where V = X(X 2 − 3Y 2 ) [14].
There the top indexes {c, v} refer to conduction and valence bands.Notice that the Γ 9 irrep appears in three pairs of basis functions, which allows for the s-p z mixing [110][111][112] Here, however, we always work on the diagonal basis (H WZ is diagonal at k = 0), which is indicated by the primes in the orbitals above.For a recent and detailed discussion on this choice of representation and the s-p z mixing, please refer to Ref. [113].
Using the basis functions from Table IV to calculate the effective 10 × 10 model using qsymm, we obtain the Hamiltonian H WZ shown in Appendix C.Here we always consider two conduction bands, which leads to this 10 × 10 generic model H WZ .However, one can also opt to work with traditional 8 × 8 models with a single conduction band.Notice, however, that for GaP the first conduction band transforms as Γ8 , while for GaN and InP the first conduction band is Γ9 .Therefore, one must be careful when selecting the appropriate 8 × 8 model for wurtzite materials.For the valence bands, one always gets Γ7 ⊕ 2 Γ9 , however, the internal ordering of these valence bands may change between materials and it can be highly sensible to the choice of density functional [28,60,114,115].The numerical coefficients c j found for GaN, GaP, InP are shown in Appendix C, and the resulting band structures are shown in Figs.4(ce).In all cases, we see that the crude model with 1000 bands (black lines) approaches a full zone description, but here we are more interested in the reduced 10 × 10 models (red lines), which present satisfactory agreement with the DFT data at low energies.

Effects of the SOC corrections on Pm,n
As introduced in Sec.II D 2, the matrix elements P m,n can be calculated with or without the PAW corrections, p SOC , that carry the SOC contributions.For most of the materials we have studied here, these corrections are marginal and the results from both cases are nearly identical.Nevertheless, we emphasize that using our patched bands.xwithin QE is faster than using the Python code to calculate P m,n via Eq.(13).
To illustrate the effects of the PAW/SOC corrections on the matrix elements P m,n , Fig. 5 compares the models for GaN and GaP with and without these corrections.For the conduction bands, we notice that the p SOC corrections significantly improve the GaN effective mass, but barely affect GaP.For the valence bands, both GaN and GaP show moderate effects of p SOC .Indeed, this shows that a precise calculation of P m,n is critical to improve the precision of the models [116].

C. Rock-salt crystals
The crystal lattice for rock-salt crystals is shown in Fig. 6(a), which is an FCC lattice with two atoms in the base, and it is described by the space group Fm 3m.The low energy band structure concentrates at the L point of the Brillouin zone shown in Fig. 6(b), which transforms as the D 3D point group after factorizing the Bloch translations.The basis functions for the first valence and conduction bands transform as where A 1g is the trivial irrep for S-like orbitals, and A 2u represent Z-like orbitals [117].Therefore, the basis functions for the L+ 6 bands are {|S, ↑⟩ , |S, ↓⟩}, and for L− 6 one gets {|Z, ↑⟩ , |Z, ↓⟩}.Here, the x, y, and z coordinates are taken along the [ 11 2], [1 10], and [111] crystallographic directions.
Here we consider two examples of rock-salt crystals: PbSe and SnTe.Their effective 4 × 4 Hamiltonian H RS under the L± 6 basis, and its numerical parameters are shown in Appendix C, and the comparison between DFT and model band structures are shown in Figs.6(c)-(d).PbSe is a narrow gap semiconductor, where the conduction band transforms as the L+ 6 irrep, and the valence band as L− 6 .In contrast, SnTe shows inverted bands, with L+ 6 below L− 6 , yielding a topological insulator phase [118,119].In both cases, the low-energy model captures the main features of the bands, including the anisotropy.

D. Other examples
To finish the set of illustrative examples, we show here the case for: (i) the monolayer MoS 2 , which is one of the most studied transition metal dichalcogenides (TMDC) [120][121][122]; (ii) the bulk bismuth selenide (Bi 2 Se 3 ), which is one of the first discovered 3D topological insulators [123,124]; and (iii) a monolayer of GaBiCl 2 , which is a large gap 2D topological insulator [125].The symmetry characteristics and basis functions for the low-energy bands of these materials mentioned above are summarized in Table V.For MoS 2 , the first valence and conduction bands are given by the single group irreps A ′ and E ′ 1 of the C 3h group [39,126], which can be represented as S-like and (X + iY )-like orbitals.For GaBiCl 2 , the valence bands are characterized by single group E irrep, and it splits into E ⊗D 1/2 = Γ4 ⊕ Γ5 ⊕ Γ6 in the spinful case, while the conduction band is given by the irrep A 1 ⊗D 1/2 = Γ6 .For Bi 2 Se 3 , a detailed derivation of the effective model can be seen in Ref. [127], which shows that the first valence and conduction bands are given by A 1g ⊗D 1/2 = Γ + 6 , and The effective Hamiltonians and their numerical coefficients for these materials can be found in the Examples folder of the code repository.Here we show only the comparison between the DFT and model band structures in Fig. 7.The MoS 2 case, as shown in Fig. 7(a), is challenging for a k • p method, since its band structure presents valleys in between high symmetry points.Consequently, the 4 bands model (red lines) captures only the nearly parabolic dispersion at the K point.However, the crude all-bands model (black lines, see Eq. ( 4)) approaches a full zone description and captures the valley along the Γ-K direction.For GaBiCl 2 , Fig. V. DISCUSSIONS Above, we have presented illustrative results of the capabilities of our code to calculate the k • p Kane and Luttinger parameters for a series of relevant materials.In all cases we see a patent agreement between the DFT (QE) data and the low-energy models near the relevant k 0 point.However, it is important to notice that here we use only PBE functionals [128], consequently it often underestimates the gap (e.g.0.5 eV instead of 1.5 eV for GaAs).Therefore, our models are limited by the quality of the DFT bands and the resulting numerical parameters might not match Kane and Luttinger's parameters for well-known materials, for which these parameters are typically chosen to match the experimental data, and not the DFT simulations.
For instance, let us consider the zincblende crystals' Kane parameter E P = 2m 0 P 2 /ℏ 2 , band gap E g and effective mass for the conduction band m * .For GaAs, the experimental values are E P ∼ 24 eV, P ∼ 0.96 eVnm, E g ∼ 1.5 eV, and m * = 0.065m 0 [58].As mentioned above, the DFT results with PBE functionals underes-timate the gap, and we get E g ∼ 0.5 eV.Moreover, the Kane parameter can be written as P = − √ 6c 5 /2, where the coefficient c 5 = −0.635eVnm is shown in Appendix C.This value yields P ∼ 0.7 eVnm and E P ∼ 16 eV.The effective mass for the conduction band can be estimated from its spinless expression [47], m 0 /m * = 1 + 2m 0 P 2 /E g ℏ 2 , which gives us m * = 0.031m 0 .While these numbers do not match well with the experimental values, we notice that if we fix the GaAs gap (scissorscut approximation), but keep our value for P , we find m * = 0.058m 0 , which is already much closer to the experimental value for the effective mass.
The number estimates shown above clearly indicate that the quality of our models is limited to the DFT simulations only.Particularly, the gap issue can be fixed if one replaces the PBE functionals with hybrid functionals, GW calculations, or other methods that improve the material gap accuracy.These are beyond the scope of this paper, but it is a possible path for future improvements of our code.
In all examples presented here, we always consider the crude all bands model from Eq. ( 4), and the optimal symmetry-adapted (few bands) model from Eq. ( 5).This raises two interesting questions: (i) how many bands are necessary for convergence?And (ii) for a large number of bands, should we get a full zone description?We discuss these questions below.

A. Convergence
The convergence threshold (how many bands are necessary) strongly depends on the material.In some cases ∼ 300 bands are sufficient, but in others, it often needs ∼ 1000 bands.We do not have a general rule to establish which materials will show a slow or fast convergence.Nevertheless, we believe it is instructive to discuss the outcomes of our convergence analysis.
Notice that the Löwdin partitioning from Eq. ( 5) has two distinct contributions.The first two terms in Eq. ( 5) are the zeroth and first-order perturbation terms.These terms do not change as we increase the number of DFT bands (provided that there are enough bands to converge the DFT calculation itself).The zeroth order term is essentially given by the DFT eigenstates, and the first order terms are given by the matrix elements ⟨m| H ′ (k) |n⟩ = 2k • P m,n between eigenstates of set A, which is the low energy sector of interest.In contrast, the third term defines the second-order corrections, which are quadratic in k (assuming a diagonal basis at k = 0).In this case, the second-order contributions depend explicitly on the sum over the remote set of bands B. These are the terms that strongly depend on the number of remote bands.
To check for convergence, we plot the values of the Hamiltonian coefficients c j associated with second-order corrections as a function of the number of remote bands.In the Examples folder in the code repository, one finds these plots for all cases presented in this paper.Here, in the top panels of Fig. 8, we select a few illustrative cases.In the bottom panels of Fig. 8 we combine the discrete derivatives of c j into a single dimensionless metric for convergence C(N ), which read as where c j (N ) refers to the coefficient calculated using N remote bands.With increasing N , the coefficients are expected to converge, consequently C(N ) → 0. The data for C(N ) is shown in blue dots on the bottom panels of Fig. 8, which is significantly noisy due to the discrete jumps on the evolution of c j with increasing N .Therefore, we also plot a moving average C(N ) (orange lines) to clearly show the convergence.For spinless graphene in Fig. 8(a), there are only two second order c j terms (neglecting terms with k z , since it is a 2D material), and we see that it reaches convergence with less than 300 remote bands.
In contrast, for MoS 2 , the convergence requires at least ∼ 500 remote bands.Interestingly, it has been recently shown that TMDC materials indeed require a large number of bands to converge the orbital angular momenta [42][43][44][45].This fact may be associated with the large number of unoccupied bands with plane-wave character that appear due to the spatial extension of the vacuum region.The GaN and GaP cases in Figs.8(c)-(d) are interesting cases, they belong to the same class of materials, but GaP reaches convergence with ∼ 200 remote bands, while GaN is not yet fully converged for ∼ 1000 remote bands.Unlike monolayer materials, the GaN compound is not described by any vacuum region, and therefore we speculate that such poor convergence may be related to details of the pseudopotential [100] and the electronegativity of Nitrogen.

B. Full zone kp
In Section II A we have presented the k•p method in its traditional form, which considers a perturbative expansion of the Bloch Hamiltonian at a reference momentum k 0 , and a small set of bands near the Fermi energy.Usually, one expects the resulting effective model to be valid only near k 0 and only for a small energy range that encloses the bands of interest.In contrast, within the full zone k•p approach [66,[129][130][131][132][133] one considers a large set of bands, such that the resulting low energy model agrees well with DFT or experimental bands over the full Brillouin zone, instead of only the vicinity of k 0 .However, to achieve this precision, one needs to apply fitting procedures to ensure that the bands match selected energy levels at various k points over the Brillouin zone.
Here, in our code, we can easily select an arbitrary number of bands to build effective models.All examples presented above show sets of bands colored in red and black, such that the red ones consider models built from a small set of bands A (from 4 to 10 bands), while the black ones consider the full set of bands from the DFT data (typically 500 or 1000 bands).This leads to an interesting question: should our all bands model match the full zone k • p models?
To answer this question, let us focus first on the graphene results from Fig. 1.There, we have seen that the QE/DFT and the model agree remarkably well at low energies near the K point, as expected.Particularly, the red line for the optimal symmetry-adapted model describes precisely the low energy regime and Dirac cone and the trigonal warping from the quadratic terms in Eq. (21).In contrast, when we consider the all-bands model (black lines), we see that the model approaches a full zone agreement with 300 bands.What if we consider more bands?Our numerical tests have shown that increasing the number of bands does improve the overall description, approaching the full zone agreement.However, this is a very slow convergence and we never really reach a true full zone agreement.This characteristic is seen in all other examples shown here.
For GaAs, Gawarecki and collaborators [133] show an excellent full zone agreement between model and DFT bands considering 30 bands.In contrast, our results presented in Fig. 3(a) for 8 (red) and 1000 (black) bands remain valid only in the vicinity of Γ.The key difference is the fitting procedure.The full zone models fit the bands over the full Brillouin zone, while in our approach we consider only the direct ab initio matrix elements of Metric C (N) Moving average Figure 8. Convergence of the second-order coefficients cj as a function of the number of remote bands for (a) spinless graphene, (b) MoS2, (c) GaN, and (d) GaP.On top (a1-d1), each panel shows the coefficients cj for different material.On panels (c1) and (d1) we omit the legends because there are 30 distinct coefficients, ranging from c22 to c51, which makes their individual identification cumbersome, and it is sufficient to visualize that all lines become nearly flat for a large number of remote bands.On the bottom (a2-d2), for each material, the evolution of the coefficients cj are combined into convergence metric set by Eq. ( 35) (blue dots).Due to the noise induced by the discrete derivative in this metric, we plot the moving average of the data as a guide for the eyes.π = p + p SOC without further manipulation.
If one needs a full zone model, we suggest using our results as the initial guess for the parameters used on a band-fitting algorithm.Moreover, since the fitted parameters must not deviate significantly from our ab initio results, our calculated values provide an important benchmark for the fitting results.Alternatively, it might be possible to develop multi-valley k • p models [68,70,134] and extract its parameters directly from DFT matrix elements without numerical fitting procedures, but this is beyond the scope of this work.

VI. CONCLUSIONS
We have implemented a numerical framework to calculate the k • p Kane and Luttinger parameters and optimal symmetry-adapted effective Hamiltonians directly from ab initio wavefunctions.The code is mostly written in Python but also contains a patch to modify the Quantum ESPRESSO code, such that its bands.xpost processing tool is used to calculate the matrix elements P m,n = ⟨m| π |n⟩, which is the central quantity in our methodology.Consequently, this first version works only with Quantum ESPRESSO.Equivalent calculations can be done in other DFT codes (e.g.VASP [5], Wien2k [6]), but it requires further developments.The code is open source and it is available at Ref. [89].
Here, we have illustrated the capabilities of our code applying it to a series of relevant and well-known materials.The resulting effective models yield band structures that match well the DFT data in the low energy sector near the k point used for the wavefunction expansion.Therefore, our code provides an ab initio approach for the k • p numerical parameters, which can be contrasted with fitting methods [60,[76][77][78]133], in which the numerical coefficients are obtained by numerically minimizing the residue difference between the DFT and model band structures over a selected range of the Brillouin zone.These fitting procedures work well in general but require careful verification if the fitted parameters are reasonable.In contrast, our ab initio approach is automatic and fully reliable.Nevertheless, fitting procedures can improve the agreement between DFT and the model band structures significantly.In this case, we suggest that our code can be used (i) to generate the initial values for the fitting parameters, and (ii) to verify if the fitted parameters show reasonable values.One should expect that fitted parameters must not deviate much from our ab initio values.
Here we do not perform a thorough comparison of our numerical parameters with experimental data.Typically, to obtain precise agreement with experimental data, one needs to fix the gap issue by using either hybrid functionals or GW calculations, which are beyond the scope of this first version of the code.Instead, here we use only PBE functionals [128] for simplicity, which is reliable enough to validate our approach.Consequently, our numerical parameters are limited by the precision of the DFT simulation, and we would not expect remarkable agreement with experimental data for most materials at this stage.Nevertheless, for novel materials, for which there is no experimental data available, our code can be used to generate reliable numerical parameters that can be improved later, either in comparison with future experiments or by extending our method to work with hybrid functionals or GW calculations.
As a final disclaimer, we would like to state that after developing the first version of the code, we have found that Ref. [86] recently proposes an equivalent approach to build k • p models from DFT, but the authors do not provide an open-source code.In any case, despite the similarities, the development of our code was done independently from their proposal.In practice, the only significant difference between the proposals is the approach to calculate the transformation matrix U (see Section II C).While the authors of Ref. [86] follow the method from [95], here we propose a different method that is more efficient for transformations involving reducible represen-tations, which is necessary when dealing with nearly degenerate bands of different irreps (e.g., spinful graphene).Additionally, after the initial submission of our paper, a new code VASP2kp [135] was released with functionalities similar to ours, but designed for VASP [5] instead of QE.
Table XI.Effective Hamiltonian for wurtzite crystals considering the 10 × 10 model with two conduction bands.

Figure 1 .
Figure 1.Graphene lattices emphasizing the Dirac cone eigenstates at the K point, where (a) |A⟩ = |(X + iY )Z⟩ and (b) |B⟩ = |(X − iY )Z⟩.Both eigenstates are composed by pz orbitals centered at the colored sites (A and B lattices) with the Bloch phase factors indicated within the circles, where τ = exp(i2π/3).(c) The first Brillouin zone, marking the path Γ − K − M used to plot the bands in (d).(d) Band structure for graphene calculated via QE/DFT (blue circles), all bands model [Eq.(4)] (black lines), and optimal symmetryadapted model [Eq.(21)] for the two bands forming the Dirac cone (red).Here, the QE/DFT simulation was performed with 300 bands.

Figure 2 .
Figure 2. The absolute value of the representation matrices of the symmetry operations for the spinful graphene example, as labeled on top of each column.The top line of matrices are defined under the ideal basis informed by the user, i.e. {|(X + iY )Z, ↑⟩, |(X − iY )Z, ↓⟩, |(X − iY )Z, ↑⟩, |(X + iY )Z, ↓⟩}, as discussed in the text.The central line shows the calculated representation matrices under the crude DFT basis from QE, which does not split into the ideal block-diagonal form due to the small SOC gap between the bands.Applying our transformation U to the crude representation from the central line, we obtain the optimal symmetryadapted basis that lead to the proper block-diagonal form of the representation matrices shown in the bottom line.

Figure 3 .
Figure 3. (a) Zincblende lattice, and (b) its first Brillouin zone (FCC).The band structure for (c) GaAs, (d) HgTe, and (e) CdTe are shown over a large energy scale on the main panels, while at the bottom of each panel, we show a zoom over the relevant low energy range.In all cases, the DFT data consider 1000 bands.

Figure 4 .
Figure 4. (a) Lattice and (b) Brillouin zone for wurtzite crystals.Band structures for (c) GaN, (d) GaP, and (e) InP show the large energy range on top, and a zoom shows the top of the valence bands at the bottom of each panel.In all cases, the DFT calculation considers 1000 bands.

Figure 5 .
Figure 5.Comparison between the DFT data and the effective models calculated with the full matrix element Pm,n including PAW/SOC corrections (red lines) and the simplified Pm,n without PAW/SOC corrections (green lines) for (a) GaN and (b) GaP.

Figure 6 .
Figure 6.(a) The rock salt lattice and (b) its Brillouin zone (FCC).Band structures for (c) PbSe and (d) SnTe.The bottom of each panel zooms into the low energy range near the Fermi level.Both DFT calculations were performed considering 500 bands.

Figure 7 .
Figure 7. Band structures for: (a) MoS2, (b) GaBiCl2, and (c) Bi2Se3 showing only the relevant low energy range.The DFT calculations were performed for 1000, 500, and 500 bands, respectively.(d) Rhombohedral lattice of Bi2Se3 and 2D hexagonal lattice of (e) MoS2 and (f) GaBiCl2, where we have omitted the vacuum region (15 Å) perpendicular to the plane formed by vectors A1 and A2.(g) 2D Brillouin zone common to MoS2 and GaBiCl2, and (h) 3D BZ of Bi2Se3.
7(c), the 6 bands model describes satisfactorily the low energy conduction and valence bands.For Bi 2 Se 3 in Fig. 7(b) the 4 bands model captures well the low-energy band structure near Γ, including the hybridization between the inverted bands.

Table I .
Graphene parameters for the Hamiltonian of Eq. (21).

Table V .
Summary of space group, irreps and basis functions for the low energy bands of MoS2, GaBiCl2, and Bi2Se3.The first column lists the materials, the second indicates the lattice space group, and the little group at the relevant k point.The third and fourth columns lists the irreps and basis functions for the low energy bands in each case.The table shows the double group irreps and the corresponding single group irreps between parenthesis.

Table X .
Effective Hamiltonian for zincblende crystals considering the 8 × 8 extended Kane model.