⇦Chinese Physics Letters, 2023, Vol. 40, No. 12, Article code 127101Express Letter VASP2KP: $k\!\cdot\! p$ Models and Landé $g$-Factors from ab initio Calculations Sheng Zhang (章盛)1,2†, Haohao Sheng (盛昊昊)1,2†, Zhi-Da Song (宋志达)3,4,5*, Chenhao Liang (梁晨昊)1,2, Yi Jiang (蒋毅)1,2, Song Sun (孙松)1,2, Quansheng Wu (吴泉生)1,2, Hongming Weng (翁红明)1,2, Zhong Fang (方忠)1,2, Xi Dai (戴希)6, and Zhijun Wang (王志俊)1,2* Affiliations 1Beijing National Laboratory for Condensed Matter Physics, and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China 2University of Chinese Academy of Sciences, Beijing 100049, China 3International Center for Quantum Materials, School of Physics, Peking University, Beijing 100871, China 4Hefei National Laboratory, Hefei 230088, China 5Collaborative Innovation Center of Quantum Matter, Beijing 100871, China 6Department of Physics, Hong Kong University of Science and Technology, Hong Kong 999077, China Received 7 November 2023; accepted manuscript online 6 December 2023; published online 13 December 2023 †These authors contributed equally to this work.
*Corresponding authors. Email: songzd@pku.edu.cn; wzj@iphy.ac.cn Citation Text: Zhang S, Sheng H H, Song Z D et al. 2023 Chin. Phys. Lett. 40 127101    Abstract The $k\!\cdot\! p$ method is significant in condensed matter physics for the compact and analytical Hamiltonian. In the presence of magnetic field, it is described by the effective Zeeman's coupling Hamiltonian with Landé $g$-factors. Here, we develop an open-source package VASP2KP (including two parts: vasp2mat and mat2kp) to compute $k\!\cdot\! p$ parameters and Landé $g$-factors directly from the wavefunctions provided by the density functional theory (DFT) as implemented in Vienna ab initio Simulation Package (VASP). First, we develop a VASP patch vasp2mat to compute matrix representations of the generalized momentum operator $\hat{\boldsymbol \pi}=\hat{\boldsymbol p}+\frac{1}{2mc^2}[\hat{{\boldsymbol s}}\times\nabla V({\boldsymbol r})]$, spin operator $\hat{\boldsymbol s}$, time reversal operator $\hat{T}$, and crystalline symmetry operators $\hat{R}$ on the DFT wavefunctions. Second, we develop a python code mat2kp to obtain the unitary transformation $U$ that rotates the degenerate DFT basis towards the standard basis, and then automatically compute the $k\!\cdot\! p$ parameters and $g$-factors. The theory and the methodology behind VASP2KP are described in detail. The matrix elements of the operators are derived comprehensively and computed correctly within the projector augmented wave method. We apply this package to some materials, e.g., Bi$_2$Se$_3$, Na$_3$Bi, Te, InAs and 1H-TMD monolayers. The obtained effective model's dispersions are in good agreement with the DFT data around the specific wave vector, and the $g$-factors are consistent with experimental data. The VASP2KP package is available at https://github.com/zjwang11/VASP2KP.
cpl-40-12-127101-fig1.png
cpl-40-12-127101-fig2.png
cpl-40-12-127101-fig3.png
cpl-40-12-127101-fig4.png
 DOI:10.1088/0256-307X/40/12/127101 © 2023 Chinese Physics Society Article Text 1. Introduction. Electronic band structures hold immense importance in the field of condensed matter physics and materials science, providing crucial insights into the behavior of electrons within crystalline materials. They reveal information about distributions of energy levels, energy gaps, and density of states, which are essential for determining material electrical conductivities, optical properties, thermal behavior, and so on. In this background, the density functional theory (DFT)[1,2] was developed, which gave birth to many first-principles calculation softwares or codes based on them, including VASP,[3,4] Quantum Espresso,[5,6] CASTEP,[7] ABINIT,[8,9] and so on.[10-13] However, DFT bands are quite complex, which contain quite a large number of energy bands that hardly affect the desired physical properties thus making the physical pictures difficult to understand. In order to construct a model that only involves a few bands with a greater impact on the physical properties of materials, the $k\!\cdot\! p$ method is proposed, which is used to construct an effective model to describe the quasiparticles near the specific wave vector in the reciprocal space.[14] The $k\!\cdot\! p$ models constructed through the theory of invariants are analytical and only contain a few important bands, thus making the physical pictures quite clear. There are some arbitrary parameters in $k\!\cdot\! p$ models, which could be determined by fitting to the corresponding experimental data or the DFT calculation data such as band structures. The $k\!\cdot\! p$ method has been successfully applied to many condensed-matter systems, including metals,[15,16] semiconductors,[17,18] topological insulators and superconductors,[19-22] spin-lasers,[23,24] nanostructured solids,[25] two-dimensional van der Waals materials,[26,27] and so on.[28-31] When a magnetic field is applied to a condensed matter system, we can use effective Zeeman's coupling Hamiltonian to describe the effects of the magnetic field. Effective Zeeman's coupling determines the split of Kramers states under magnetic field, leading to outcomes like the Pauli paramagnetism of the metals and the van Vleck paramagnetism of insulators. The important parameters in Zeeman's coupling are Landé $g$-factors. Landé $g$-factors of materials have been widely studied in quantum wires,[32-35] quantum dots,[36,37] semiconductors nanostructures,[38-42] topological materials,[43,44] and so on.[45-47] Moreover, the theory of invariants can also be used to construct Zeeman's coupling Hamiltonian conveniently. However, there is no available code to compute Landé $g$-factors and $k\!\cdot\! p$ parameters directly from DFT wavefunctions. First of all, Gao et al. developed the program IRVSP to determine the irreducible representations (irreps) of bands at any $\boldsymbol k$-point in VASP calculations.[48] Then, Jiang et al. developed the python package kdotp-generator to generate the $k\!\cdot\! p$ effective Hamiltonian for the given irreps automatically.[49] Furthermore, Song et al. derived effective Zeeman's coupling and $g$-factors in DFT calculations.[50] Thus, it is straightforward for us to generate the VASP2KP code to construct the effective models and to compute the model parameters from VASP wavefunctions directly. Although many similar codes are developed subsequently for Quantum Espresso, such as IR2PW,[51] IrRep,[52] and DFT2kp,[53] their functions do not go beyond the above-mentioned codes. The DFT2kp can not generate Zeeman's coupling Hamiltonians. Additionally, the matrix elements of symmetry operators are not obtained properly, because the projector augmented wave (PAW) corrections are neglected in DFT2kp. In this work, we have developed an open-source package VASP2KP, which can generate the effective $k\!\cdot\! p$ model and Zeeman's coupling, and obtain the values of the $k\!\cdot \! p $ parameters and $g$-factors. This package contains two parts: a VASP patch vasp2mat and a post-processing python code mat2kp. First, we use vasp2mat to generate matrix representations for generalized momentum $\hat{\boldsymbol \pi}=\hat{\boldsymbol p}+\frac{1}{2mc^2}(\hat{\boldsymbol s}\times\nabla V(\boldsymbol r))$, spin $\hat{\boldsymbol s}$, time reversal $\hat{T} $, and crystalline symmetry operators $\hat{R}$ in DFT calculations. The matrix elements of the space group operators are derived in detail and computed correctly in the PAW wavefunctions. Second, mat2kp can obtain the unitary transformation $U $ from the degenerate wavefunctions to the $k\!\cdot\! p$ standard basis, and then compute these parameters automatically. The obtained effective masses and $g$-factors are important and comparable with experimental observations. In the following, we first introduce the theoretical foundations behind VASP2KP. Then, the main algorithm steps are described in detail. The results obtained by VASP2KP are shown and analyzed for some typical materials. Finally, we present some discussions. 2. Theory and Methodology. Now, we first introduce the $k\!\cdot \! p $ method. The method to obtain Zeeman's coupling is presented. Then the theory of invariants is reviewed and the invariant $k\!\cdot\! p$ Hamiltonian as well as Zeeman's coupling are obtained. We propose a general routine to get the unitary transformation that changes degenerate DFT basis to the $k \cdot p$ standard basis. Lastly, the method to calculate $k \cdot p$ parameters and $g$-factors is presented. 2.1 $k$ $\cdot$ $p$ effective Hamiltonian. When we only care about wave vector ${\boldsymbol K} $ around a specific wave vector ${\boldsymbol k}_0 $ in the Brillouin zone, it can be well described by using the $k\!\cdot \! p $ effective model. Suppose that $\psi_{n{\boldsymbol K}}({\boldsymbol r}) $ is the wavefunction of $n$-th band which satisfies the Schrödinger equation with spin-orbit coupling (SOC), i.e., $\hat{H}_{\scriptscriptstyle{\rm B}}\psi_{n{\boldsymbol K}}({\boldsymbol r})=\epsilon_{n}({\boldsymbol K})\psi_{n{\boldsymbol K}}({\boldsymbol r}) $, where \begin{align} \hat{H}_{\scriptscriptstyle{\rm B}}=\frac{\hat{p}^2}{2\,m}+V({\boldsymbol r})+\frac{1}{2m^2c^2}[\hat{{\boldsymbol s}}\times\nabla V({\boldsymbol r})]\cdot\hat{{\boldsymbol p}} \tag {1} \end{align} is the Bloch Hamiltonian operator with SOC, and $\hat{\boldsymbol p} $ is the momentum operator, $V({\boldsymbol r}) $ is the potential in crystal, $\hat{{\boldsymbol s}}$ is the spin momentum operator, $m$ is the electron mass, and $c$ is the light speed in vacuum. To expand the Hamiltonian at $\boldsymbol k_0$, introducing a transformation $\phi_{n{\boldsymbol k}}({\boldsymbol r})=\psi_{n{\boldsymbol K}}({\boldsymbol r})e^{-i\boldsymbol k\cdot\boldsymbol r} $ with ${\boldsymbol k}={\boldsymbol K}-{\boldsymbol k}_0$ being the deviation from ${\boldsymbol k}_0$, the Schrödinger equation that $\phi_{\boldsymbol k}({\boldsymbol r})$ obeys is \begin{align} \Big[\hat{H}_{\scriptscriptstyle{\rm B}}+\frac{\hbar^2k^2}{2\,m}+\hat{H}'({\boldsymbol k})\Big]\phi_{n{\boldsymbol k}}({\boldsymbol r})=\epsilon_n({\boldsymbol K})\phi_{n{\boldsymbol k}}({\boldsymbol r}), \tag {2} \end{align} where $\hat{H}'({\boldsymbol k})=\frac{\hbar}{m}{\boldsymbol k}\cdot{\boldsymbol \pi}$ is the first-order term, and $\hat{\boldsymbol \pi}=\hat{\boldsymbol p}+\frac{1}{2mc^2}[\hat{{\boldsymbol s}}\times\nabla V({\boldsymbol r})]$ is the generalized momentum operator with SOC. We can take $\hat{H}^{kp}=\hat{H}_{\scriptscriptstyle{\rm B}}+\frac{\hbar^2k^2}{2\,m}+\hat{H}'({\boldsymbol k})$ as the equivalent Hamiltonian of $\hat{H}_{\scriptscriptstyle{\rm B}}$ since the eigenvalues of them are all equal. Suppose that we have got the eigenenergies $\{\epsilon_n({\boldsymbol k}_0)\}$ and the eigenstates $\{|n({\boldsymbol k}_0)\rangle\} $ of $\hat{H}_{\scriptscriptstyle{\rm B}}$. Then we can obtain the generalized momentum elements by ${\boldsymbol \pi}_{mn}=\langle m({\boldsymbol k}_0)|{\hat{\boldsymbol \pi}}|{n({\boldsymbol k}_0)}\rangle$. Thus, the matrix elements of $\hat{H}^{kp} $ can be obtained by \begin{align} H^{kp{\textrm{-}}{\rm all}}_{mn}({\boldsymbol k})=\Big[\epsilon_n({\boldsymbol k}_0)+\frac{\hbar^2k^2}{2\,m}\Big]\delta_{mn}+\frac{\hbar}{m}{\boldsymbol \pi}_{mn}\cdot{\boldsymbol k}. \tag {3} \end{align} Usually, we aim at several low-energy bands (the set of which is denoted as $\mathcal{A}$). The set of other bands is denoted as $\mathcal{B}$. Then we can fold down the $H^{kp{\textrm{-}}{\rm all}}_{mn}({\boldsymbol k})$ Hamiltonian into a subspace $\mathcal{A}$ via Löwdin partitioning. After two-order Löwdin partitioning, the effective $k\!\cdot \! p $ Hamiltonian is transformed as (see Section A in the Supplementary Material for details) \begin{align} H^{kp}_{\alpha\beta}({\boldsymbol k})=\,&\Big[\epsilon_\alpha({\boldsymbol k}_0)+\frac{\hbar^2k^2}{2\,m}\Big]\delta_{\alpha\beta}+\frac{\hbar}{m}{\boldsymbol \pi}_{\alpha\beta}\cdot{\boldsymbol k}\notag\\ &+\frac{\hbar^2}{2m^2}\sum_{l\in \mathcal{B}}\sum_{ij}\pi^{i}_{\alpha l}\pi^{j}_{l\beta}k^ik^j\notag\\ &\times \Big[\frac{1}{\epsilon_\alpha({\boldsymbol k}_0)-\epsilon_l({\boldsymbol k}_0)}+\frac{1}{\epsilon_\beta({\boldsymbol k}_0)-\epsilon_l({\boldsymbol k}_0)}\Big] \tag {4} \end{align} with $i,j\in \{x,y,z\} $. Hereafter we use $\alpha, \beta\in\mathcal{A}$, $l\in\mathcal{B}$ and $m, n\in \mathcal{A}\cup \mathcal{B}$. Moreover, the $k \cdot p$ Hamiltonian of the third order can be constructed similarly in Section A of the Supplementary Material. 2.2 Zeeman's Coupling. When a magnetic field is applied to the system, the momenta $\hbar k^i$ should be replaced by $(-i\hbar\partial^i+eA^i)$ according to Peierls substitution, where ${\boldsymbol A}$ is the vector potential and $e$ (positively valued) is the elementary charge. Therefore, $\hbar^2k^ik^j $ in the summation will be replaced by the sum of the gauge dependent term $\frac{1}{2}\{-i\hbar\partial^i+eA^i,-i\hbar\partial^j+eA^j\}$ and gauge independent term $\frac{1}{2}[-i\hbar\partial^i+eA^i,-i\hbar\partial^j+eA^j]=-\frac{i\hbar e}{2}\sum_k\epsilon^{ijk}B_k$, where ${\boldsymbol B}$ is the magnetic induction intensity. In this case, the total Hamiltonian can be written by the summation of the $k \cdot p$ effective Hamiltonian $H^{kp}_{\alpha\beta} $ and Zeeman's coupling $H^{Z}_{\alpha\beta}$, which are gauge dependent and gauge invariant, respectively. The $k \cdot p$ effective Hamiltonian $H^{kp}_{\alpha\beta}$ can be expressed as \begin{align} H^{kp}_{\alpha\beta}=\,&\epsilon_\alpha({\boldsymbol k}_0)\delta_{\alpha\beta}+\frac{\hbar}{m}{\boldsymbol \pi}_{\alpha\beta}\cdot\Big(-i\nabla+\frac{e}{\hbar}{\boldsymbol A}\Big)\notag\\ &+\sum_{ij}M_{\alpha\beta}^{ij}\Big(-i\partial^i+\frac{e}{\hbar}A^i\Big) \Big(-i\partial^j+\frac{e}{\hbar}A^j\Big), \tag {5} \end{align} where \begin{align} M_{\alpha\beta}^{ij}=\,&\frac{\hbar^2}{2\,m}\delta_{\alpha\beta}\delta_{ij}+\frac{\hbar^2}{4m^2}\sum_{l\in\mathcal{B}}(\pi^{i}_{\alpha l}\pi^{j}_{l\beta}+\pi^{j}_{\alpha l}\pi^{i}_{l\beta})\notag\\ &\times \Big[\frac{1}{\epsilon_\alpha({\boldsymbol k}_0)-\epsilon_l({\boldsymbol k}_0)}+\frac{1}{\epsilon_\beta({\boldsymbol k}_0)-\epsilon_l({\boldsymbol k}_0)}\Big]. \tag {6} \end{align} Zeeman's coupling term can be expressed by \begin{align} H^{Z}_{\alpha\beta}=\frac{\mu_{\scriptscriptstyle{\rm B}}}{\hbar}({\boldsymbol L}_{\alpha\beta}+2{\boldsymbol s}_{\alpha\beta})\cdot{\boldsymbol B}, \tag {7} \end{align} where \begin{align} L_{\alpha\beta}^{k}=\,& - \frac{i\hbar}{2\,m}\sum_{l\in\mathcal{B}}\sum_{ij}\epsilon^{ijk}\pi^{i}_{\alpha l}\pi^{j}_{l\beta}\notag\\ &\times\Big[\frac{1}{\epsilon_\alpha({\boldsymbol k}_0)-\epsilon_l({\boldsymbol k}_0)}+\frac{1}{\epsilon_\beta({\boldsymbol k}_0)-\epsilon_l({\boldsymbol k}_0)}\Big] \tag {8} \end{align} can be considered as the orbital contribution, ${\boldsymbol s}_{\alpha\beta}=\langle\alpha({\boldsymbol k}_0)|{\hat{\boldsymbol s}}|{\beta({\boldsymbol k}_0)}\rangle$ are the spin elements, and $\mu_{\scriptscriptstyle{\rm B}}$ is the Bohr magneton. The detailed derivation is shown in Section B of the Supplementary Material. 2.3 Theory of Invariants. Suppose that the little group at ${\boldsymbol k}_0 $ is $L$, and $D(R) $ is a representation (rep) for $\forall R\in L$, whose dimension is $N$. There are a total of $N^2$ independent Hermitian matrices ${\boldsymbol X}={X_i}$ $(i=1,\,2,\,\ldots,\,N^2)$ in $N$ dimensions that constitute an $N$-dimensional matrix space. It is simple to get a new rep $D^{(X)}$ by \begin{align} D(R)X_nD(R)^{-1}=\sum_{m}X_mD^{(X)}_{mn}(R) \tag {9} \end{align} for $\forall R\in L $. Also, we can construct polynomial space of order $p\in\mathbb{N} $, whose basis can be chosen as ${\boldsymbol g}({\boldsymbol k})=\{k_x^{i}k_y^{j}k_z^{l}|i+j+l\leq p,i,j,l\in \mathbb{N}\} $. It is simple to get a new rep $D^{(g)}$ by \begin{align} \hat{R}g_n(\boldsymbol k)=g_n(R^{-1}{\boldsymbol k})=\sum_m g_m({\boldsymbol k})D^{(g)}_{mn}(R). \tag {10} \end{align} Spatial inversion $\hat{P}$ and time reversal $\hat{T}$ may also be the elements of the little group $L$, which transform the wave vector as $\boldsymbol k \xrightarrow{\hat{P}} -\boldsymbol k$ and $\boldsymbol k \xrightarrow{\hat{T}} -\boldsymbol k$, respectively. Decomposing $D^{(X)} $ and $D^{(g)} $ into irreps, we have \begin{align} \tag {11} &D^{(X)}(R)\simeq\oplus_\lambda a_\lambda D^{\lambda}(R),\\ \tag {12} &D^{(g)}(R)\simeq\oplus_\lambda b_\lambda D^{\lambda}(R), \end{align} where $D^{\lambda}$ represents irreps, and $a_\lambda$ or $b_\lambda$ denotes the multiplicities of $D^{\lambda}$. Note that the irreps decomposed from the $D^{(X)}$ and $D^{(g)}$ are not exactly the same. It is simple to get the matrix basis $X^{\lambda,\mu}$ and polynomial basis $\boldsymbol g^{\lambda,\nu}(\boldsymbol k)$ of $D^{\lambda}$-irrep. An irrep may correspond to more than one sets of matrix basis or polynomial basis, thus enabling us to add an indicator $\mu/\nu$ to distinguish them. After obtaining matrix basis and polynomial basis, the $k \cdot p$ Hamiltonian can be written as \begin{align} H^{kp}({\boldsymbol k})=\sum_{\lambda}\sum_{\mu=1}^{a_\lambda}\sum_{\nu=1}^{b_\lambda}c_{\lambda\mu\nu}\boldsymbol X^{\lambda,\mu}\cdot\boldsymbol g^{*\lambda,\nu}({\boldsymbol k}), \tag {13} \end{align} where $c_{\lambda\mu\nu} $ are undetermined $k\!\cdot\! p$ parameters, which must be real.[49] It can be proved that the $p$-th order $k\!\cdot\! p$ Hamiltonian expressed by Eq. (13) satisfies the relation \begin{align} H^{kp}(\hat{R}{\boldsymbol k})=D(R)H^{kp}({\boldsymbol k})D(R)^{-1} \tag {14} \end{align} for $\forall R\in L$.[49] Moreover, in the presence of magnetic field, Zeeman's coupling can also be constructed in the same way. Unlike the wave vector $\boldsymbol k$, the magnetic field $\boldsymbol B$ is a pseudovector, thus making it satisfy the relations $\boldsymbol B \xrightarrow{\hat{P}} \boldsymbol B $ and $\boldsymbol B \xrightarrow{\hat{T}} -\boldsymbol B$ under spatial inversion and time reversal, respectively. Effective Zeeman's coupling satisfies the relation \begin{align} H^Z(\hat{R}\boldsymbol B) = D(R)H^Z(\boldsymbol B)D(R)^{-1}. \tag {15} \end{align} By the way, given the matrix representations of the generators of the little group $L$, the python package kdotp-generator constructs the standard $k\!\cdot\! p$ Hamiltonian or Zeeman's coupling with undetermined $k\!\cdot\! p$ parameters as Eq. (13). The standard matrices $D^{{\rm std}}(R) $ are given on the Bilbao Crystalline Server (BCS) for the assigned irreps by IRVSP. 2.4 Unitary Transformation. In DFT calculations, when the irrep is $n$-fold ($n\geq 2$) and the eigenstates are degenerate, there is $U(n)$ ambiguity in the DFT eigenstates. We propose a general routine to get the unitary transformation that changes degenerate DFT wavefunctions to the $k\!\cdot\! p$ standard basis. We can obtain the eigenenergies $\{\epsilon_n({\boldsymbol k}_0)\}$ and eigenstates $\{|n({\boldsymbol k}_0)\rangle\}$ in VASP. The matrix representation of $\hat{R}$ under the VASP basis set would be calculated by mat2kp, which is denoted as $D^{{\rm num}}(R) $. However, $D^{{\rm num}}(R) $ is usually not in the standard form; they are related by a unitary transformation, i.e., \begin{align} D^{{\rm std}}(R) = U^†D^{{\rm num}}(R)U \tag {16} \end{align} for $\forall R \in L$, where $U$ is a unitary matrix. To find $U$, it is obvious that only the generators of the $\boldsymbol k_0$-little group $L$ should be taken into account. They are space group operators $\{R_t|\boldsymbol v_t\}$, consisting of rotational part and translational part. The translational part is usually expressed by a phase factor. Equation (16) can be rewritten as \begin{align} D^{{\rm num}}(R)U-UD^{{\rm std}}(R)=\boldsymbol 0, \tag {17} \end{align} where $\boldsymbol 0$ is a zero matrix. The real part and the imaginary part of each elements of $U$ are independent variables, which are denoted as $U_{r11}, U_{r12}, \ldots U_{rnn} $ and $U_{i11}, U_{i12}, \ldots U_{inn} $, respectively. From Eq. (17), it is clear to find that all these variables satisfy a linear equation set. Using a column vector ${\boldsymbol u}=(U_{r11},\, U_{r12},\, \ldots,\,U_{rnn},\,U_{i11},\,U_{i12},\, \ldots,\,U_{inn})^{\scriptscriptstyle{\rm T}} $, Eq. (17) can be rewritten as \begin{align} Q{\boldsymbol u}=\boldsymbol 0, \tag {18} \end{align} where $Q $ is a real coefficient matrix defined by Eq. (17). The details are given in Section C of the Supplementary Material. Therefore, all vectors in the null space of the matrix $Q$ is the solution to Eq. (16). Performing singular value decomposition on matrix $Q$, we can obtain $Q=V_1\varSigma V_2^{\scriptscriptstyle{\rm T}}$ thus $QV_2=V_1\varSigma$, where $V_1$ and $V_2$ are orthogonal matrices and $\varSigma$ is a diagonal matrix. It can be proved that the column vectors in $V_1$ correspond to the singular value 0 forming the basis of the null space of the matrix $Q$, which are denoted as $\{{\boldsymbol u_1},{\boldsymbol u_2}, \ldots,{\boldsymbol u_m}\}$. After reshaping these vectors, the basis of the solution space of $U$ of Eq. (16) are denoted as $\{U_1,U_2,\ldots,U_n\}$. Any solution can be written as \begin{align} U=\sum_{\alpha=1}^{m}\lambda_\alpha U_\alpha , \tag {19} \end{align} where $\lambda_\alpha$ must be real. Thus, one has to get one set of the parameters $\{\lambda_\alpha\}$ to satisfy $U^†U=\mathbb{I}$, where $\mathbb{I}$ is an identity matrix. However, this relation is nonlinear, thus making it hard to solve by treating it as an equation directly. In mat2kp, the unitary $U$ is obtained via the sequential least squares programming to find optimal parameter set $\{\lambda_\alpha\}$, which can minimize the error $\varepsilon=\sum_{ij}|(U^†U-\mathbb{I})_{ij}|$. 2.5 Calculations of Parameters. Applying Eqs. (4–8) in DFT calculations, one can get the numerical Hamiltonian $H^{\rm num}$, one have to solve the equation $H^{{\rm std}}=U^{-1}H^{{\rm num}}U$ to get all the parameters in $H^{\rm std}$. By the way, from Eq. (13), it is easy to find that the equations are all linear equations. However, there are more equations than parameters. The coefficients of each linear equation are not accurate because of numerical errors generated from VASP calculations, making it impossible to solve by the traditional Gaussian elimination method. We can write all equations in a matrix form \begin{align} A\boldsymbol{x}=\boldsymbol{b}, \tag {20} \end{align} where $A $ is a constant matrix, and $\boldsymbol{x}$ is the column vector comprised of all undetermined real parameters. To eliminate the constraint of the real $\boldsymbol{x}$, suppose that $A=A_r+iA_i $ and $\boldsymbol{b}=\boldsymbol{b}_r+i\boldsymbol{b}_i$, where the subscript $r$ represents the real part and the subscript $i$ represents the imaginary part. Then Eq. (20) can be transformed as \begin{align} \begin{pmatrix} A_r\\A_i \end{pmatrix}\boldsymbol{x}=\begin{pmatrix} \boldsymbol{b}_r\\ \boldsymbol{b}_i \end{pmatrix}. \tag {21} \end{align} In this case, all the matrices become real (the real constraint is automatically satisfied). Therefore, the linear least square method can be used to calculate the parameters and to give the error as well.
cpl-40-12-127101-fig1.png
Fig. 1. The workflow of VASP2KP to compute the $k \cdot p$ parameters and Landé $g$-factors directly from the VASP wavefunctions.
3. Capability of VASP2KP. The VASP2KP package contains two parts: the VASP patch vasp2mat and the post-processing python code mat2kp. The workflow is presented in Fig. 1. The methodology of obtaining the matrix elements of the generalized momentum, spin, time reversal and crystalline rotational operators in VASP is introduced in the following. The main algorithm steps of VASP2KP are also described. Lastly, main steps to construct $k \cdot p$ models and Zeeman's coupling via VASP2KP are present. 3.1 vasp2mat: to Compute Matrix Elements from DFT Wavefunctions. The matrix elements of $\hat{\boldsymbol \pi}$, $\hat{\boldsymbol s}$, $\hat{T}$, and $\hat{R}$ are required in constructing the $k \cdot p$ Hamiltonian and Zeeman's coupling under the DFT wavefunctions as shown in Section 2.1. The computation method of $\hat{\boldsymbol \pi}$ matrix elements was first introduced in Ref. [50]. Here we generalize it to nonlocal operators. In VASP, the PAW method is used, which is a combination and generalization of the linear augmented plane wave method. Therefore, we derive these matrix elements under the PAW wavefunctions, which are implemented in vasp2mat to compute the matrix elements. The relation of all electron wavefunction ($|n(\boldsymbol k_0)\rangle$) and pseudo wavefunction $(|\widetilde{n}(\boldsymbol k_0)\rangle)$ in PAW is \begin{align} |n(\boldsymbol k_0)\rangle =\mathcal{T}|\widetilde{n}(\boldsymbol k_0)\rangle, \tag {22} \end{align} where the linear transformation $\mathcal{T}$ can be expressed as[54,55] \begin{align} \mathcal{T}=1+\sum_{a\mu\zeta}\Big(|\phi^a_{\mu}\zeta\rangle-|\widetilde{\phi}^a_{\mu}\zeta\rangle\Big)\langle\widetilde{p}^a_{\mu}\zeta|. \tag {23} \end{align} The $|\phi^a_{\mu}\zeta\rangle$ is the direct product of the real space all electron partial wavefunction $\phi^a_{\mu}(\boldsymbol r)$ at the $a$-th atom with a spinor wavefunction $|\zeta\rangle(\zeta=\uparrow {\rm or}\downarrow)$, $|\widetilde{\phi}^a_{\mu}\zeta\rangle$ is the direct product of the real space pseudo partial wavefunction $\widetilde{\phi}^a_{\mu}(\boldsymbol r)$ with a spinor wavefunction $|\zeta\rangle$, and $|\widetilde{p}^a_{\mu}\zeta\rangle$ is a projector wavefunction comprised of the direct product of a real space projector wavefunction $\widetilde{p}^a_{\mu}({\boldsymbol r})$ and a spinor wavefunction $|\zeta\rangle$. Here, $\phi^{a}_{\mu}(\boldsymbol r)$'s are obtained by all-electron calculation for the reference atom. The pseudo partial wavefunctions $\widetilde{\phi}^a_{\mu}(\boldsymbol r)$ are identical to $\phi^{a}_{\mu}(\boldsymbol r)$ outside the augmentation sphere of the corresponding atom and are much softer than $\phi^{a}_{\mu}(\boldsymbol r)$ inside the augmentation sphere; $\widetilde{\phi}^{a}_{\mu}(\boldsymbol r)$'s provide a complete basis for the pseudo wavefunction $|\widetilde{n}(\boldsymbol{k}_0)\rangle$ inside the augmentation sphere. The projector wavefunctions are defined in such a way, i.e., $\langle \widetilde{p}^a_{\mu} | \widetilde{\phi}^{a}_{\mu'} \rangle=\delta_{\mu\mu'}$, that $\langle \widetilde{p}_{\mu}^a\zeta | \widetilde{n}(\boldsymbol{k}_0)\rangle$ gives the expanding coefficients of $|\widetilde{n}(\boldsymbol{k}_0)\rangle$ on $\widetilde{\phi}^a_{\mu}(\boldsymbol r)$. Usually the projector wavefunctions are chosen to be zero outside the core radius. Therefore, outside the augmentation sphere the third and second terms in $\mathcal{T}$ cancel each other exactly, and inside the augmentation sphere the first and third terms cancel each other exactly. $\mathcal{T}$ leaves $| \widetilde{n}(\boldsymbol{k}_0)\rangle$ unchanged outside the augmentation sphere and maps it to the all-electron wavefunction inside the augmentation sphere. Both ${\phi}^a_{\mu}(\boldsymbol r) $ and $\widetilde{\phi}^a_{\mu}(\boldsymbol r)$ are stored as a radial part times an angular part, which can be expressed as \begin{align} &\phi^{a}_{\mu}(\boldsymbol r)=Y^{m_\mu}_{l_\mu}(\widehat{\boldsymbol r-{\boldsymbol \tau}_a})R^a_{\mu}(|\boldsymbol r-{\boldsymbol \tau}_a|),\notag\\ &\widetilde{\phi}^{a}_{\mu}(\boldsymbol r)=Y^{m_\mu}_{l_\mu}(\widehat{\boldsymbol r-{\boldsymbol \tau}_a})\widetilde{R}^a_{\mu}(|\boldsymbol r-{\boldsymbol \tau}_a|), \tag {24} \end{align} where ${\boldsymbol \tau}_a$ is the site of the $a$-th atom, $Y^{m}_l$ are sphere harmonics, and $R^{a}_\mu$ are real functions. The hatted vectors represent the unit vector along the corresponding directions. According to Eq. (23), the PAW matrix form of a local operator $\hat{F}$ can be expressed by \begin{align} F_{mn}=\,&\langle m(\boldsymbol k_0)|{\hat{F}}|{n(\boldsymbol k_0)}\rangle=\langle{\widetilde{m}(\boldsymbol k_0)}|{\mathcal{T}^†\hat{F}\mathcal{T}}|{\widetilde{n}(\boldsymbol k_0)}\rangle\notag\\ =\,&\langle{\widetilde{m}(\boldsymbol k_0)}|{\hat{F}}|{\widetilde{n}(\boldsymbol k_0)}\rangle\notag\\ &+\sum_{a\mu\nu}\sum_{\zeta\zeta'}\langle{\widetilde{m}(\boldsymbol k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle F^{a}_{\mu\zeta,\nu\zeta'}\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\boldsymbol k_0)}\rangle, \tag {25} \end{align} where $F^{a}_{\mu\zeta,\nu\zeta'}$ is the projection matrix of the operator $\hat{F}$ in the $a$-th atom's augmentation sphere, which is defined as \begin{align} F^{a}_{\mu\zeta,\nu\zeta'}=\langle{\phi^a_{\mu}\zeta}|{\hat{F}}|{\phi^{a}_{\nu}\zeta'}\rangle -\langle{\widetilde{\phi}^a_{\mu}\zeta}|{\hat{F}}|{\widetilde{\phi}^{a}_{\nu}\zeta'}\rangle. \tag {26} \end{align} The first term in $F^{a}_{\mu\zeta,\mu'\zeta'}$ gives the contribution from the all-electron wavefunction in the augmentation sphere, and the second term cancels the contribution from the pseudo wavefunction in the augmentation sphere that is counted in the first term of Eq. (25). In PAW, the plane wave is used to span the pseudo wavefunction $\widetilde{\psi}_{n \boldsymbol k_0}(\boldsymbol r,\,\zeta)=\langle{\boldsymbol r,\zeta}|{\widetilde{n}(\boldsymbol k_0)}\rangle$, which is expressed by \begin{align} \widetilde{\psi}_{n\boldsymbol k_0}(\boldsymbol r,\, \zeta)=\sum_{\boldsymbol G}c^{n\boldsymbol k_0}_{\zeta\boldsymbol G}e^{i({\boldsymbol k_0+\boldsymbol G})\cdot{\boldsymbol r}}, \tag {27} \end{align} where $c^{n\boldsymbol k_0}_{\zeta\boldsymbol G} $ are the plane wave coefficients. The projection coefficients $\langle{\widetilde{p}^a_{\mu}\zeta}|{\widetilde{n}(\boldsymbol k_0)}\rangle$ in the second term of Eq. (25) have already been calculated in VASP. Therefore, to calculate the matrix of the local operator $\hat{F}$, we need to calculate $\langle{\widetilde{m}(\boldsymbol k_0)}|{\hat{F}}|{\widetilde{n}(\boldsymbol k_0)}\rangle$ and $F^{a}_{\mu\zeta,\nu\zeta'}$ in Eq. (25). Generalized Momentum Matrix. The matrix of the generalized momentum $\hat{\boldsymbol \pi}$ can be expressed by substituting $\hat{\boldsymbol \pi}$ into $\hat{F}$ in Eq. (25) as follows: \begin{align} &\boldsymbol \pi_{mn}=\langle{\widetilde{m}(\boldsymbol k_0)}|{\hat{\boldsymbol p}}|{\widetilde{n}(\boldsymbol k_0)}\rangle+\frac{1}{2mc^2}\langle{\widetilde{m}(\boldsymbol k_0)}|{\hat{\boldsymbol s}\times \nabla V}|{\widetilde{n}(\boldsymbol k_0)}\rangle\notag\\ +&\sum_{a\mu\nu}\sum_{\zeta\zeta'}\langle{\widetilde{m}(\boldsymbol k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\boldsymbol p^{a}_{\mu\zeta,\nu\zeta'}\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\boldsymbol k_0)}\rangle\notag\\ +&\frac{1}{2mc^2}\sum_{a\mu\nu\zeta\zeta'}\langle{\widetilde{m}(\boldsymbol k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\langle{\phi^a_{\mu}\zeta}|{\hat{\boldsymbol s}\times\nabla V}|{\phi^{a}_{\nu}\zeta'}\rangle\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\boldsymbol k_0)}\rangle\notag\\ -&\frac{1}{2mc^2}\sum_{a\mu\nu\zeta\zeta'}\langle{\widetilde{m}(\boldsymbol k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\langle{\widetilde{\phi}^a_{\mu}\zeta}|{\hat{\boldsymbol s}\times\nabla V}|{\widetilde{\phi}^{a}_{\nu}\zeta'}\rangle\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\boldsymbol k_0)}\rangle. \tag {28} \end{align} Since the SOC effect is considered only within the augmentation spheres in VASP, where $|\widetilde{n}(\boldsymbol k_0)\rangle =\sum_{a\mu\zeta}|\widetilde{\phi}^a_{\mu}\zeta\rangle \langle{\widetilde{p}^a_\mu\zeta}|{\widetilde{n}(\boldsymbol k_0)}\rangle$, we can make the second term and the fifth term in Eq. (28) cancel out. Therefore, the generalized momentum matrix can be simplified to \begin{align} \boldsymbol \pi_{mn}=&\langle{\widetilde{m}(\boldsymbol k_0)}|{\hat{\boldsymbol p}}|{\widetilde{n}(\boldsymbol k_0)}\rangle\notag\\ &+\sum_{a\mu\nu}\sum_{\zeta\zeta'}\langle{\widetilde{m}(\boldsymbol k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle\boldsymbol \pi'^{a}_{\mu\zeta,\nu\zeta'}\langle{\widetilde{p}^a_{\nu}\zeta'}|{\widetilde{n}(\boldsymbol k_0)}\rangle, \tag {29} \end{align} where \begin{align} \boldsymbol \pi'^{a}_{\mu\zeta,\nu\zeta'}=\,&\delta_{\zeta\zeta'}\Big(\langle{\phi^a_\mu}|{\hat{{\boldsymbol p}}}|{\phi^a_\nu}\rangle-\langle{\widetilde{\phi}^a_\mu}|{\hat{{\boldsymbol p}}}|{\widetilde{\phi}^a_\nu}\rangle\Big)\notag\\ &+\frac{\hbar^2}{4mc^2}\boldsymbol \sigma_{\zeta\zeta'}\times\langle{\phi^a_\mu}|{\nabla V}|{\phi^a_\nu}\rangle. \tag {30} \end{align} The first term of $\boldsymbol \pi_{mn} $ in Eq. (29) can be calculated by \begin{align} \langle{\widetilde{m}(\boldsymbol k_0)}|{\hat{\boldsymbol p}}|{\widetilde{n}(\boldsymbol k_0)}\rangle=\sum_{\boldsymbol G}\hbar(\boldsymbol k_0+\boldsymbol G)c^{*m\boldsymbol k_0}_{\zeta\boldsymbol G}c^{n\boldsymbol k_0}_{\zeta\boldsymbol G}, \tag {31} \end{align} where $\hat{\boldsymbol p}=-i\hbar\nabla$. The integrals in $\boldsymbol \pi'^a_{\mu\zeta,\nu\zeta'}$ can be calculated by separating the radial part and angular part, which are expressed as \begin{align} \langle{\phi^a_\mu}|{\hat{{\boldsymbol p}}}|{\phi^a_\nu}\rangle=\,&-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\nabla Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2R^{a*}_\mu R^a_\nu\notag\\ &-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\frac{\boldsymbol r}{r}Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2R^{a*}_{\mu}\partial_rR^{a}_\nu,\notag\\ \langle{\widetilde{\phi}^a_\mu}|{\hat{{\boldsymbol p}}}|{\widetilde{\phi}^a_\nu}\rangle=\,&-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\nabla Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2\widetilde{R}^{a*}_\mu \widetilde{R}^a_\nu\notag\\ &-i\hbar\int\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\frac{\boldsymbol r}{r}Y^{m_\nu}_{l_\nu}\int\mathrm{d}rr^2\widetilde{R}^{a*}_{\mu}\partial_r\widetilde{R}^{a}_\nu,\notag\\ \langle{\phi^a_\mu}|{\nabla V}|{\phi^a_\nu}\rangle\approx\,&\!\int\!\mathrm{d}\Omega Y^{m_\mu*}_{l_\mu}\frac{\boldsymbol r}{r}Y^{m_\nu}_{l_\mu}\int\mathrm{d}rr^2R^{a*}_{\mu}\partial_rVR^{a}_\nu. \tag {32} \end{align} Finally, by substituting Eqs. (30)-(32) into Eq. (29), the matrices of the generalized momentum $\hat{\boldsymbol \pi} $ can be obtained. Spin Matrices. By substituting the local operator $\hat{\boldsymbol s} $ for $\hat{F} $ in Eq. (25), the corresponding matrix elements are given explicitly. The first term in Eq. (25) can be calculated by \begin{align} \langle{\widetilde{\alpha}(\boldsymbol k_0)}|{\hat{\boldsymbol s}}|{\widetilde{\beta}(\boldsymbol k_0)}\rangle=\frac{\hbar}{2}\sum_{\boldsymbol G\boldsymbol G'\zeta\zeta'}\boldsymbol \sigma_{\zeta\zeta'}c^{*\alpha\boldsymbol k_0}_{\zeta\boldsymbol G}c^{\beta\boldsymbol k_0}_{\zeta'\boldsymbol G'}\delta_{\boldsymbol G,\boldsymbol G'}, \tag {33} \end{align} where $\hat{\boldsymbol s}=\frac{\hbar}{2}\hat{\boldsymbol \sigma}$ and the projection matrix can be obtained by \begin{align} \!{\boldsymbol s}^{a}_{\mu\zeta,\nu\zeta'}=&\frac{\hbar}{2}\delta_{l_\mu l_\nu}\delta_{m_\mu m_\nu}\boldsymbol \sigma_{\zeta\zeta'}\int\mathrm{d}rr^2(R^{a*}_\mu R^{a}_\nu-\widetilde{R}^{a*}_\mu \widetilde{R}^{a}_\nu). \tag {34} \end{align} Finally, by substituting Eqs. (33) and (34) into Eq. (25), the matrices of the spin $\hat{\boldsymbol s}$ can be obtained in the PAW wavefunctions. Space Group Operator Matrices. The general (magnetic) space group operator (SGO) is expressed by $F=\{R|\boldsymbol v\}$ or $F=\{TR|\boldsymbol v\}$, which is usually a nonlocal (NL) operator. One should notice that the SGO commutes with $\mathcal{T}$ in Eq. (23). Thus the matrix element can be written as \begin{align} F^{NL}_{mn}=\,& \langle{m(\boldsymbol k_0)}|{\hat{F}}|{n(\boldsymbol k_0)}\rangle=\langle{\widetilde{m}(\boldsymbol k_0)}|{\mathcal{T}^† \mathcal{T}}|{\hat{F} \cdot \widetilde{n}(\boldsymbol k_0)}\rangle \notag\\ =\,& \langle{\widetilde{m}(\boldsymbol k_0)}|\hat{F} \cdot \widetilde{n}(\boldsymbol k_0)\rangle \notag\\ &+\sum_{a\mu\nu\zeta}\langle{\widetilde{m}(\boldsymbol k_0)}|{\widetilde{p}^a_{\mu}\zeta}\rangle C^{a}_{\mu,\nu}\langle{\widetilde{p}^{a}_{\nu}\zeta}|{\hat{F} \cdot \widetilde{n}(\boldsymbol k_0)}\rangle, \tag {35} \end{align} where \begin{align} &C^a_{\mu,\nu}=\delta_{l_\mu l_\nu}\delta_{m_\mu m_\nu}\int\mathrm{d}rr^2(R^{a*}_\mu R^{a}_\nu-\widetilde{R}^{a*}_\mu \widetilde{R}^{a}_\nu),\notag\\ &\langle{\boldsymbol r,\zeta}|{\{R|\boldsymbol v\}}|{\widetilde{\beta}(\boldsymbol k_0)}\rangle=\sum_Gc^{\beta\boldsymbol k_0}_{\zeta\boldsymbol G} e^{i({\boldsymbol k_0}+{\boldsymbol G})\cdot R^{-1}({\boldsymbol r}-{\boldsymbol v})},\notag\\ &\langle{\boldsymbol r,\zeta}|{\{TR|\boldsymbol v\}}|{\widetilde{\beta}(\boldsymbol k_0)}\rangle=\eta_\zeta\sum_Gc^{\beta\boldsymbol k_0*}_{\overline{\zeta}\boldsymbol G} e^{-i({\boldsymbol k_0}+{\boldsymbol G})\cdot R^{-1}({\boldsymbol r}-{\boldsymbol v})}. \tag {36} \end{align} Here, $\eta_\uparrow=1$, $\eta_\downarrow=-1$, and $\bar\zeta$ indicates the opposite of spin $\zeta$. In our convention, the translation $\boldsymbol L$ is expressed by a phase factor of $e^{-i \boldsymbol k_0\cdot \boldsymbol L}$. More functions of the patch vasp2mat are presented in Section D of the Supplementary Material. 3.2 The Python Code mat2kp. In this subsection, we give a brief introduction of the main algorithm steps of mat2kp. This code needs the inputs of $\varepsilon_n(\boldsymbol k_0)$, $\boldsymbol {\pi}_{mn}$, $\boldsymbol {s}_{\alpha\beta}$, $D^{{\rm num}}_{\alpha\beta}(R)$, and $D^{{\rm std}}_{\alpha\beta}(R)$, which correspond to the EIGENVAL, MAT_Pi.m, MAT_sig.m, MAT_R.m, and mat2kp.in files, respectively. The main steps of mat2kp are as follows:
  1. Following Section 2.4, calculate the unitary transformation matrix $U $ satisfying $D^{{\rm std}}(R) = U^† D^{{\rm num}}(R)U $ for $\forall R\in L$.
  2. By downfolding processing and using Eqs. (5)-(8), get the numerical $k\!\cdot \! p $ Hamiltonian and Zeeman's coupling, and compute the coefficients numerically under DFT basis, i.e., $H^{kp-{\rm num}}$ and $H^{Z-{\rm num}}$.
  3. Through theory of invariants in Section 2.3, import the package kdotp-generator and generate the standard $k\!\cdot \! p $ Hamiltonian and Zeeman's coupling with a set of undetermined parameters, $H^{kp-{\rm std}} $ and $H^{Z-{\rm std}} $, respectively.
  4. Obtain the values of $k\!\cdot \! p $ parameters and Landé $g $-factors by solving $H^{kp-{\rm std}}=U^{-1}H^{kp-{\rm num}}U $ and $H^{Z-{\rm std}}=U^{-1}H^{Z-{\rm num}} U$ in Section 2.5.
 3.3 General Steps to Get the $k\!\cdot\! p$ Model and Parameters Automatically. The general workflow is given in Fig. 1. We will take Bi$_2$Se$_3$ as an example for illustration.
  1. Run VASP to output the eigenstates $|n(\boldsymbol k_0)\rangle$ $(|\widetilde{n}(\boldsymbol k_0)\rangle$ in WAVECAR) and eigenvalues $\varepsilon_n(\boldsymbol {k}_0)$ (EIGENVAL) at $\boldsymbol k_0$ point.
  2. Run IRVSP to get the irreps of the aimed low-energy bands (set $\mathcal{A} $), and then obtain the standard matrix representations [$D^{{\rm std}}(R) $] of the generators of $\boldsymbol k_0$-little group on the Bilbao Crystalline Server (BCS). They are given in “mat2kp.in” (the input file of mat2kp).
    ######## mat2kp.in - Bi2Se3 ########
    Symmetry = {
    'C3z' : {
    'rotation_matrix':
    Matrix([[-Rational(1,2), -sqrt(3)/2,0],[sqrt(3)/2, -Rational(1,2), 0],[0, 0, 1]]),
    'repr_matrix':
    Matrix([[Rational(1,2)-I*sqrt(3)/2,0,0,0],[0,Rational(1,2)+I*sqrt(3)/2,0,0],
    [0,0,Rational(1,2)-I*sqrt(3)/2,0],[0,0,0,Rational(1,2)+I*sqrt(3)/2]]),
    'repr_has_cc': False},
    'C2x' : {
    'rotation_matrix': Matrix([[1, 0, 0],[0, -1, 0],[0, 0, -1]]),
    'repr_matrix':
    Matrix([[0,-Rational(1,2)-I*sqrt(3)/2,0,0],[Rational(1,2)-I*sqrt(3)/2,0,0,0],
    [0,0,0,-Rational(1,2)-I*sqrt(3)/2],[0,0,Rational(1,2)-I*sqrt(3)/2,0]]),
    'repr_has_cc': False},
    'P' : {
    'rotation_matrix': Matrix([[-1,0,0],[0, -1, 0],[0, 0, -1]]),
    'repr_matrix':
    Matrix([[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]]),
    'repr_has_cc': False},
    'T' : {
    'rotation_matrix': eye(3),# Identity Matrix
    'repr_matrix': Matrix([[0,1,0,0],[-1,0,0,0],[0,0,0,-1],[0,0,1,0]]),
    'repr_has_cc': True}
    }
    # optional parameters
    vaspMAT = '../Bi2Se3/GMmat' # the path: to read eigenvalues, Pi, s, and R matrices in this folder.
    order = 2     # Order of the kp model : 2 (default) or 3.
    print_flag = 2 # Where to output results: 1 (screen) or 2 (files, default).
    kpmodel = 1   # Whether to compute Hkp: 0 or 1 (default).
    gfactor = 1   # Whether to compute HZ : 0 or 1 (default).
    log = 1     # Whether to output log files: 0 or 1 (default).

  3. Run vasp2mat to generate $\boldsymbol \pi_{mn}$ (vma$t=11$; vmat_name='Pi'), $\boldsymbol s_{\alpha\beta}$ (vma$t=10$; vmat_name='sig'), and generators' [$D(R)^{{\rm num}}_{\alpha\beta}$; vma$t=12$] matrices directly from the VASP wavefunctions (WAVECAR) by the following settings in “INCAR.mat” files, respectively. These numerical matrix representations are output in MAT_Pi.m, MAT_sig.m, MAT_R.m files. They are in the ‘vaspMAT’ folder (given in mat2kp.in). The vmat_name of the generators should be the same as those given in “mat2kp.in”.
    ######## INCAR.mat - Pi mat. ########
    &vmat_para
    ! soc————————
    cfactor=1.0
    socfactor=1.0
    nosoc_inH = .false.
    ! operator——————
    vmat = 11
    vmat_name = 'Pi'
    vmat_k = 1
    bstart=1, bend=400
    print_only_diagnal = .false.
    /

    ######## INCAR.mat - sigma mat. ########
    &vmat_para
    ! soc————————
    cfactor=1.0
    socfactor=1.0
    nosoc_inH = .false.
    ! operator——————
    vmat = 10
    vmat_name = 'sig'
    vmat_k = 1
    bstart=47, bend=50
    print_only_diagnal = .false.
    /

    ######## INCAR.mat - C3z / C2x / P / T mat. - Bi2Se3 ########
    &vmat_para
    ! soc————————
    cfactor=1.0
    socfactor=1.0
    nosoc_inH = .false.
    ! operator——————
    vmat = 12
    vmat_name = 'C3z'  /  'C2x'  /  'P'   /  'T'
    vmat_k = 1
    bstart=47, bend=50
    print_only_diagnal = .false.
    ! rotation——————
    rot_n(:) = 0 0 1   /  1 0 0  /  0 0 1  /  0 0 1
    rot_alpha = 120    /  180     /  0     /  0
    rot_det = 1      /  1     /  -1    /  1
    rot_tau(:) = 0 0 0
    rot_spin2pi = .false.
    time_rev = .false.  /  .false. /  .false. /  .true.
    /

  4. Run the post-processing python code mat2kp with the input files, to construct the standard $k\!\cdot \! p $ Hamiltonian and Zeeman's coupling, and compute the values of $k\!\cdot\! p$ parameters and $g$-factors. The outputs are “kp-parameters.out” and “g-factors.out” files, as pasted below.
    ######## kp-parameters.out ########
    kp Hamiltonian
    ========== Result of kp Hamiltonian ==========
    Matrix([[a1 + a2 + c1*(kx**2 + ky**2) + c2*(kx**2 + ky**2) + c3*kz**2 + c4*kz**2, 0, -I*b2*kz, -b1*(kx*(sqrt(3) + 3*I) + ky*(3 - sqrt(3)*I))/3], [0, a1 + a2 + c1*(kx**2 + ky**2) + c2*(kx**2 + ky**2) + c3*kz**2 + c4*kz**2, b1*(kx*(sqrt(3) - 3*I) + ky*(3 + sqrt(3)*I))/3, I*b2*kz], [I*b2*kz, b1*(kx*(sqrt(3) + 3*I) + ky*(3 - sqrt(3)*I))/3, a1 - a2 + c1*(kx**2 + ky**2) - c2*(kx**2 + ky**2) + c3*kz**2 - c4*kz**2, 0], [-b1*(kx*(sqrt(3) - 3*I) + ky*(3 + sqrt(3)*I))/3, -I*b2*kz, 0, a1 - a2 + c1*(kx**2 + ky**2) - c2*(kx**2 + ky**2) + c3*kz**2 - c4*kz**2]])
    Parameters:
    a1 = 4.8898 ;
    a2 = -0.2244 ;
    b1 = -3.238 ;
    b2 = 2.5562 ;
    c1 = 19.5842 ;
    c2 = 44.4746 ;
    c3 = 1.8117 ;
    c4 = 9.5034 ;
    Error of the linear least square method: 3.93e-06
    Sum of absolute values of numerical zero elements: 6.47e-02

    ######## g-factors.out ########
    Zeeman's coupling
    ========== Result of Zeeman's coupling ==========
    mu_B/2*Matrix([[Bz*g3 + Bz*g4, g1*(Bx*(1 - sqrt(3)*I/3) + By*(-sqrt(3)/3 - I)) + g2*(Bx*(1 - sqrt(3)*I/3) + By*(-sqrt(3)/3 - I)), 0, 0], [g1*(Bx*(1 + sqrt(3)*I/3) + By*(-sqrt(3)/3 + I)) + g2*(Bx*(1 + sqrt(3)*I/3) + By*(-sqrt(3)/3 + I)), -Bz*g3 - Bz*g4, 0, 0], [0, 0, Bz*g3 - Bz*g4, g1*(Bx*(1 - sqrt(3)*I/3) + By*(-sqrt(3)/3 - I)) + g2*(Bx*(-1 + sqrt(3)*I/3) + By*(sqrt(3)/3 + I))], [0, 0, g1*(Bx*(1 + sqrt(3)*I/3) + By*(-sqrt(3)/3 + I)) + g2*(Bx*(-1 - sqrt(3)*I/3) + By*(sqrt(3)/3 - I)), -Bz*g3 + Bz*g4]])
    Parameters:
    g1 = -0.3244 ;
    g2 = 5.761 ;
    g3 = -7.8904 ;
    g4 = -13.0138 ;
    Error of the linear least square method: 6.11e-08
    Sum of absolute values of numerical zero elements: 4.12e-03

cpl-40-12-127101-fig2.png
Fig. 2. Crystal structure and electronic band structure (vasp) of Bi$_2$Se$_3$ (a), Na$_3$Bi (c), and Te (e). Model dispersions of Bi$_2$Se$_3$ (b), Na$_3$Bi (d), and Te (f) [$kp$; Eq. (5) and $kp$-all; Eq. (3)].
4. Applications in Materials. We apply this package to some typical materials, i.e., Bi$_2$Se$_3$, Na$_3$Bi, Te, InAs, and 1H-TMD, to construct the effective models and to compute all the parameters. 4.1 Four-Band Model at $\varGamma$ in Bi$_2$Se$_3$. As we know, the topological property of Bi$_2$Se$_3$ is due to the band inversion at $\varGamma$. By using IRVSP, the lowest conduction band belongs to $\overline{{\rm GM}}9$ (twofold degenerate), while the highest valence band belongs to $\overline{{\rm GM}}8$ (twofold degenerate), as depicted in Fig. 2(a). Based on these states in the ascending order, the low-energy effective Hamiltonian at $\varGamma$ is constructed automatically. The generators of the $\varGamma$-little group are $\{C_{3z}|0,0,0\}$, $\{C_{2x}|0,0,0\}$, $\{P|0,0,0\} $, and $\{T|0,0,0\} $. Their standard matrix representations are given in Table 1, which are needed by the code mat2kp to construct effective models. To the second order, the $k \cdot p$ Hamiltonian and Zeeman's coupling are given in Eq. (37), where $\mu_{\scriptscriptstyle{\rm B}}$ is the Bohr magneton ($\sim$ $0.05788$ meV/Tesla). After their numerical matrix representations are computed by vasp2mat directly from the VASP wavefunctions, the parameters of $k\!\cdot\! p$ Hamiltonian and $g$-factors of Zeeman's coupling are computed, as shown in Table 2. The four-band model's dispersions are plotted in Fig. 2(b). They fit well with the VASP bands in the vicinity of $\varGamma$. Moreover, the dispersions of the all-band $k\!\cdot\! p$ model without downfolding in Eq. (3) are also plotted for comparison (labeled by $kp$-all). \begin{align} &H^{{\rm eff}}(\boldsymbol k,\,\boldsymbol B)= H^{kp}+H^Z,\notag\\ &H^{kp}=\begin{pmatrix} D_1&0&-ib_2k_z&-\frac{3i+\sqrt{3}}{3}b_1k_-\\ &D_1&\frac{\sqrt{3}-3i}{3}b_1k_+&ib_2k_z\\ &&D_2&0\\†&&&D_2 \end{pmatrix},\notag\\ &H^Z=\frac{\mu_{\scriptscriptstyle{\rm B}}}{2}\begin{pmatrix} h_{1}^{+}B_z&h_{2}^{+}B_+&0&0\\&-h_{1}^{+}B_z&0&0\\&&h_{1}^{-}B_z&h_{2}^{-}B_-\\†&&&-h_{1}^{-}B_z \end{pmatrix},\\ &D_1=a_1+a_2+(c_1+c_2)k_+k_-+(c_3+c_4)k_z^2,\notag\\ &D_2=a_1-a_2+(c_1-c_2)k_+k_-+(c_3-c_4)k_z^2,\notag\\ &h_{1}^{\pm}=g_3\pm g_4,\ h_{2}^{\pm}=\frac{3-\sqrt{3}i}{3}(g_1\pm g_2),\notag\\ &k_\pm=k_x\pm ik_y, \ B_\pm=B_x\pm iB_y.\notag \tag {37} \end{align}
Table 1. The matrix representations of ${\rm\overline{GM}}8$ and $\rm{\overline{GM}}9$ irreps at $\varGamma$ are given on BCS server for the generators.
$\overline{{\rm GM}}8 $ $\overline{{\rm GM}}9 $
$\{C_{3z}|0,0,0\} $ $\begin{pmatrix} e^{-\frac{\pi i}{3}} & 0\\0 & e^{\frac{\pi i}{3}} \end{pmatrix}$ $\begin{pmatrix} e^{-\frac{\pi i}{3}} & 0\\0 & e^{\frac{\pi i}{3}} \end{pmatrix}$
$\{C_{2x}|0,0,0\} $ $\begin{pmatrix} 0 & e^{-\frac{2\pi i}{3}}\\ e^{-\frac{\pi i}{3}} & 0 \end{pmatrix}$ $\begin{pmatrix} 0 & e^{-\frac{2\pi i}{3}}\\ e^{-\frac{\pi i}{3}} & 0 \end{pmatrix}$
$\{P|0,0,0\} $ $\begin{pmatrix} 1 & 0\\0 & 1 \end{pmatrix}$ $\begin{pmatrix} -1 & 0\\0 & -1 \end{pmatrix}$
$\{T|0,0,0\} $ $\begin{pmatrix} 0 & 1\\-1 & 0 \end{pmatrix}\mathcal{K}$ $\begin{pmatrix} 0 & -1\\1 & 0 \end{pmatrix}\mathcal{K}$
Table 2. The computed values of parameters $\{a_i, b_i, c_i, g_i\} $ for Bi$_2$Se$_3$ in Eq. (37), obtained from the VASP calculations directly.
$a_i$ (eV) $b_i$ (eV$\cdot$Å) $c_i$ (eV$\cdot$Å$^2$) $g_i$
$a_1=4.89$ $b_1=3.24$ $c_1=19.58$ $g_1=-0.32$
$a_2=-0.22$ $b_2=-2.56$ $c_2=44.47$ $g_2=5.76$
$c_3=1.81$ $g_3=-7.90$
$c_4=9.50$ $g_4=-13.01$
4.2 Four-Band Model at the Dirac Point in Na$_3$Bi. The Dirac semimetal Na$_3$Bi has a Dirac point (DP: $\boldsymbol k_D$) along $\varGamma$–$A$ in Fig. 2(c). It is formed by the crossing of the $\rm\overline{DT}7$-irrep and $\rm\overline{DT}8$-irrep bands. Thus, we construct the effective Hamiltonian at $\boldsymbol k_D$. The standard matrix representations are presented in Table E1 in the Supplementary Material. The $H^{kp}$ and $H^Z$ are obtained in Eq. (38) with the computed parameters and $g$-factors in Table 3. The four-band $k\!\cdot\! p$ model's dispersions are fitting well with the VASP bands, as shown in Fig. 2(d). \begin{align*} H^{kp}=&\begin{pmatrix} D_1&0&\xi_-c_2k_-^2&\varTheta_{14}\\ &D_1&-i\xi_-k_+(b_1+c_4k_z)&\xi_+c_2k_+^2\\ &&D_2&0\\ †&&&D_2 \end{pmatrix}, \end{align*} \begin{align*} &D_1=a_1+a_2+(b_2+b_3)k_z+d_1^{+}k_+k_-+d_2^+k_z^2,\notag\\ &D_2=a_1-a_2+(b_2-b_3)k_z+d_1^-k_+k_-+d_2^-k_z^2,\notag\\ &d_1^\pm=c_1\pm c_3,\notag\\ &d_2^\pm=c_5\pm c_6,\notag\\ &\xi_\pm=1\pm\sqrt{3}i,\notag\\ &\varTheta_{14}=-i\xi_+k_-(b_1+c_4k_z), \end{align*} \begin{align} H^Z=&\frac{\mu_{\scriptscriptstyle{\rm B}}}{2}\begin{pmatrix} h_{1}^{+}B_z&0&0&\varUpsilon_{14}\\ &-h_{1}^{+}B_z&\frac{3+\varUpsilon_{23}\sqrt{3}i}{3}g_2B_-&0\\ &&h_{1}^{-}B_z&\varUpsilon_{34}\\ †&&&-h_{1}^{-}B_z \end{pmatrix},\\ h_{1}^{\pm}=\,&g_3\pm g_4,~\varUpsilon_{14}=\frac{3\!-\!\sqrt{3}i}{3}g_2B_-, ~\varUpsilon_{34}=\frac{6\!+\!2\sqrt{3}i}{3}g_1B_+. \notag \tag {38} \end{align}
Table 3. The computed values of parameters $\{a_i, b_i, c_i, g_i\} $ in Eq. (38) for Na$_3$Bi, obtained from the VASP calculations directly.
$a_i$ (eV) $b_i$ (eV$\cdot$Å) $c_i$ (eV$\cdot$Å$^2$) $g_i$
$a_1=2.22$ $b_1=0.99$ $c_1=7.39$ $g_1=-3.78$
$a_2=0.00$ $b_2=1.48$ $c_2=-4.02$ $g_2=-5.24$
$b_3=-1.69$ $c_3=-0.76$ $g_3=2.74$
$c_4=7.24$ $g_4=8.58$
$c_5=3.17$
$c_6=-4.43$
4.3 Four-Band Model at H in Te. Element tellurium is a narrow-gap semiconductor. The direct gap is at $H$. The low-energy bands at $H$ are the $\rm\overline{H}4$-irrep and $\rm\overline{H}5$-irrep valence bands and $\rm\overline{H}6$-irrep (doubly degenerate) conduction bands. The four-band effective model is constructed accordingly. The standard matrix representations of the generators of the $H$-little group are given in Table E2 of the Supplementary Material. The $k\!\cdot\! p$ model at $H$ is expressed as \begin{align*} H^{kp}_{11}=\,&a_{1} + 2 a_{2} + a_{3} + (c_{1}+2c_3+c_8)k_-k_+ \notag\\ &+ (c_{14} +2 c_{15}+ c_{16})k_{z}^{2},\notag\\ H^{kp}_{12}=\,&2 (b_{6} - i b_{7}) k_{z},\notag\\ H^{kp}_{13}=\,&-\frac{\sqrt{3}}{6}[(2i\xi_+b_2+\sqrt{3}i\xi_+ b_3+\sqrt{3}\xi_+ b_4+4 ib_5)k_+\notag\\ &+(2i\xi_+ c_{10}+\sqrt{3}i\xi_+c_{11}-\sqrt{3}\xi_+c_{12}-4 ic_{13})k_+k_z,\notag\\ &+(2i\xi_+ c_4+\sqrt{3}i\xi_+c_5+\sqrt{3}\xi_+c_6+4 ic_7)k_-^2],\notag\\ H^{kp}_{14}=\,&\frac{\sqrt{3}}{3}[(2 b_2+\sqrt{3} b_3-\sqrt{3} i b_4+\xi_- b_5)k_-\notag\\ &-(2 c_{10}-\sqrt{3}c_{11}+\sqrt{3}i c_{12}-\xi_-c_{13})k_-k_z\notag\\ &+(2c_{4}+\sqrt{3}c_{5}-\sqrt{3}ic_{6}+\xi_-c_{7})k_+^2],\notag\\ H^{kp}_{22}=\,&a_{1} - 2 a_{2} + a_{3} + (c_{1}-2c_3+c_8)k_-k_+ \notag\\ &+(c_{14}- 2 c_{15}+ c_{16})k_{z}^{2},\\ H^{kp}_{23}=\,&\frac{\sqrt{3}}{3}[(-2b_{2}+\sqrt{3}b_{3}+\sqrt{3}i b_{4}-\xi_+b_{5})k_+\notag\\ &+(2c_{4}-\sqrt{3}c_{5}-\sqrt{3}ic_{6}+\xi_+c_{7})k_-^2\\ &+(2 c_{10}+\sqrt{3}c_{11}+\sqrt{3}ic_{12}-\xi_+c_{13})k_+k_z], \end{align*} \begin{align} H^{kp}_{24}=\,&\frac{\sqrt{3}}{6}[(2i\xi_-b_{2}-\sqrt{3} i\xi_-b_{3}+\sqrt{3}\xi_-b_{4}+4ib_{5})k_-\notag\\ &+\!(2i\xi_-c_{10}\!+\!\sqrt{3}i\xi_-c_{11}\!-\!\sqrt{3}\xi_-c_{12}\!-\!4ic_{13})k_-k_z\notag\\ &+(2i\xi_-c_{4}-\sqrt{3}i\zeta_-c_{5}-\sqrt{3}\xi_-c_{6} +4ic_{7})k_+^2],\notag\\ H^{kp}_{33}=\,&H^{kp}_{44}=a_1-a_3+2 b_8k_z+(c_1-c_8)k_-k_+\notag\\ &+(c_{14}-c_{16})k_z^2,\notag\\ H^{kp}_{34}=\,&\frac{2\sqrt{3}}{3}i\xi_-b_{1}k_+\!+\!\frac{2\sqrt{3}}{3} i\xi_-c_2k_-^2\!+\!2\xi_-c_{9}k_+k_z,\\ \xi_\pm=\,&1\pm\sqrt{3}\, i. \notag \tag {39} \end{align} Zeeman's coupling is expressed by \begin{align} H^Z=&\frac{\mu_{\scriptscriptstyle{\rm B}}}{2}\begin{pmatrix} 0&(g_6-ig_7)B_z&h_1B_+&h_2B_-\\ &0&h_3B_+&h_4B_+\\ &&g_8B_z&\frac{6+2\sqrt{3}i}{3}g_1B_+\\ †&&&-g_8B_z \end{pmatrix}, \notag\\ h_1=&\frac{3-\sqrt{3}i}{3}g_{2}+\frac{\sqrt{3}-i}{2}g_{3}-\frac{\sqrt{3}i+1}{2}g_{4}-\frac{2\sqrt{3}}{3}i g_{5},\notag\\ h_2=&\frac{2\sqrt{3}}{3}g_2+g_3-ig_4+\frac{\sqrt{3}-3i}{3}g_5,\notag\\ h_3=&-\frac{2\sqrt{3}}{3}g_2+g_3+ig_4-\frac{\sqrt{3}+3i}{3}g_5,\notag\\ h_4=&\frac{3+\sqrt{3}i}{3}g_2-\frac{\sqrt{3}+i}{2}g_3-\frac{\sqrt{3}i-1}{2}g_4+\frac{2\sqrt{3}}{3}ig_5. \tag {40} \end{align} The computed parameters $\{a_i, b_i, c_i, g_i\}$ are presented in Table 4. The four-band $k\!\cdot\! p$ model's dispersions agree well with the VASP bands in Fig. 2(f). 4.4 Two-Band Model at $\varGamma$ in InAs. Here, we compute the band structure of the wurtzite (WZ) InAs semiconductor. Since InAs is usually of $n$-type, we consider the $\overline{\rm GM}8$ conduction bands (doubly degenerate) for the electron doped samples. The matrix representations of generators are presented in Table E3 of the Supplementary Material. The two-conduction-band effective model is constructed in Eq. (41), and the $k\!\cdot\! p$ parameters and $g$-factors are computed as listed in Table 5. In Fig. 3(f), the splitting of the conduction bands at field 10 T is 6.7 meV, indicating an effective $g$-factor of $11.6$. It is consistent with the experimental value in the bulk material.[40,56]
Table 4. The computed values of parameters $\{a_i, b_i, c_i, g_i\} $ in Eqs. (39) and (40) for Te, obtained from the VASP calculations directly.
$a$ (eV) $b$ (eV$\cdot$Å) $c$ (eV$\cdot$Å$^2$) $g$
$a_1=5.72$ $b_1=-0.16$ $c_1=6.13$ $g_1=1.52$
$a_2=-0.03$ $b_2=0.84$ $c_2=-0.03$ $g_2=-3.12$
$a_3=-0.13$ $b_3=-1.32$ $c_3=2.51$ $g_3=-2.64$
$b_4=-1.30$ $c_4=4.16$ $g_4=0.18$
$b_5=-0.55$ $c_5=1.63$ $g_5=4.08$
$b_6=0.65$ $c_6=-7.48$ $g_6=-0.08$
$b_7=-0.98$ $c_7=-1.72$ $g_7=-8.92$
$b_8=-0.29$ $c_8=-5.96$ $g_8=-10.84$
$c_9=1.41$
$c_{10}=-2.56$
$c_{11}=5.52$
$c_{12}=-6.40$
$c_{13}=-1.04$
$c_{14}=8.47$
$c_{15}=-0.09$
$c_{16}=-46.46$
cpl-40-12-127101-fig3.png
Fig. 3. Crystal structure, electronic structure, and $k\!\cdot\! p$ band dispersions of InAs: (a) crystal structure, (b) electronic band structure (VASP), (c) model dispersions [$kp$; Eq. (5) and $kp$-all; Eq. (3)], and (d)–(f) the Zeeman effect with magnetic field along $\hat{x}$, $\hat{y}$, and $\hat{z}$, respectively. The splitting due to the magnetic field is depicted explicitly.
In addition, another effective model for the six valence bands is constructed in Section F of the Supplementary Material. Surprisingly, the highest valence bands show a splitting of 19.8 meV under 10 T magnetic field in Fig. 3(f), indicating a remarkable effective $g$-factor, being three times of that of the conduction bands. \begin{align} H^{kp}=&a_1+\begin{pmatrix} c_1k_+k_-+c_2k_z^2&\left(1+\sqrt{3}i\right)b_1k_-\\ †&c_1k_+k_-+c_2k_z^2 \end{pmatrix},\notag\\ H^{Z}=&\frac{\mu_{\scriptscriptstyle{\rm B}}}{2}\begin{pmatrix} g_2B_z&\frac{3-\sqrt{3}i}{3}g_1B_-\\ †&-g_2B_z \end{pmatrix}. \tag {41} \end{align}
Table 5. The computed values of parameters $\{a_i, b_i, c_i, g_i\} $ in Eq. (41) for wurtzite InAs are obtained from the VASP calculations directly.
$a$ (eV) $b$ (eV$\cdot$Å) $c$ (eV$\cdot$Å$^2$) $g$
$a_1=4.37$ $b_1=-0.20$ $c_1=124.41$ $g_1=-7.66$
$c_2=123.40$ $g_2=-11.60$
4.5 Two-Band Model in 1H-TMD Monolayers. In 1H-phase transition metal chalcogenide (TMD) monolayers, their direct gaps are at $K$. The two valence bands at $K$ belong to the $\rm\overline{K}8$-irrep and $\rm\overline{K}11$-irrep, respectively. The standard matrix representations of the generators of the $K$-little group are given in Table E4 of the Supplementary Material. The two-band effective models are constructed in Eq. (42) (to the second order), which are plotted in Fig. 4. The computed $k\!\cdot\! p$ parameters and $g$-factors are listed in Table 6. \begin{align} H^{kp}=&a_1+\begin{pmatrix} a_2+(c_1+c_2)k_-k_+&0\\ 0&-a_2+(c_1-c_2)k_-k_+ \end{pmatrix},\notag\\ H^Z=&\frac{\mu_{\scriptscriptstyle{\rm B}}}{2}B_z\begin{pmatrix} g_1+g_2&0\\ 0&g_1-g_2 \end{pmatrix}. \tag {42} \end{align}
Table 6. The computed values of parameters $\{a_i, b_i, c_i, g_i\} $ in Eq. (42) for 1H-TMD monolayers, obtained from the VASP calculations.
MoTe$_2$ MoSe$_2$ WTe$_2$ WSe$_2$
$a_1$ (eV) $-1.92$ $-2.93$ $-1.79$ $-2.75$
$a_2$ (eV) $-0.11$ $-0.09$ $-0.24$ $-0.23$
$c_1$ (eV$\cdot$Å$^2$) $-3.83$ $-3.89$ $-8.06$ $-6.48$
$c_2$ (eV$\cdot$Å$^2$) $0.60$ $0.47$ $2.75$ $1.73$
$g_1$ $6.82$ $6.92$ $8.94$ $8.16$
$g_2$ $-2.34$ $-2.26$ $-3.40$ $-2.88$
5. Discussion. In this work, we develop an open-source package VASP2KP[57] to construct the $k \cdot p$ Hamiltonian and Zeeman's coupling, and to compute the $k \cdot p$ parameters and Landé $g$-factors from the VASP calculations directly. By applying this package in many typical materials, we get $k\!\cdot \! p $ effective models, whose band dispersions are in good agreement with the VASP data around the specific wave vector in the Brillouin zone. As the orbital contribution between the low-energy bands is usually remarkable (due to the small energy difference), the computed $g$-factor is not accurate in the multiple-band model. Thus, we recompute the effective $g$-factor for the two conduction bands only in Bi$_2$Se$_3$ and Te. The recomputed values are mainly consistent with the previous experimental data listed in Table 7. We reveal that the conduction bands of Te have the larger $g$-factor than the valence bands, and the 1st valence band has the larger $g$ factor than the 2nd valence band in 1H-TMD monolayers. Our results indicate that VASP2KP has the capabilities to handle more materials.
Table 7. The computed effective $g$-factors for the two conduction/valence bands only in Bi$_2$Se$_3$, Te, WZ-InAs, and 1H-TMD. The conduction bands of Te have the larger $g$-factor than the valence bands. The 1st valence (v$_1$) band has the larger $g$-factor than the 2nd valence band (v$_2$) in 1H-TMD.
Band $|g_z|$ $|g_\perp|$ Experiment
Bi$_2$Se$_3$ c $20.52$ $14.9$ $32$,$23$[58]
Te c $6.40$ $3.04$
WZ-InAs c $11.60$ $7.66$ $13$[59]
MoTe$_2$ v$_1$/v$_2$ $9.17/4.49$
MoSe$_2$ v$_1$/v$_2$ $9.19/4.67$
WTe$_2$ v$_1$/v$_2$ $12.34/5.55$
WSe$_2$ v$_1$/v$_2$ $11.04/5.27$ $12.2$[60]
WSe$_2$ c$_1$/c$_2$ $1.36/6.95$ $1.72/7.68$[60]
cpl-40-12-127101-fig4.png
Fig. 4. Electronic structures and $k\!\cdot\! p$ band dispersions of 1H-TMD monolayers: (a) MoTe$_2$, (b) MoSe$_2$, (c) WTe$_2$, (d) WSe$_2$.
The obtained models provide the effective masses and $g$-factors, which are important physical quantities of the materials. The minor discrepancy from the experimental data is due to the non-precise band gap in the DFT calculations. To improve these parameters, one can use our code in the hybrid functionals or $GW$ calculations. For the new synthesized or predicted materials, for which there is no experimental data available, our code can be used to predict reliable parameters. In conclusion, VASP2KP would be widely used in the Materials Science. Acknowledgements. This work was supported by the National Key R&D Program of Chain (Grant No. 2022YFA1403800), the National Natural Science Foundation of China (Grant Nos. 11974395, 12188101, 11925408, 12274436, and 11921004), the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB33000000), and the Center for Materials Genome. Zhi-Da Song was supported by the Innovation Program for Quantum Science and Technology (Grant No. 2021ZD0302403), the National Natural Science Foundation of China (Grant No. 12274005), and the National Key Research and Development Program of China (Grant No. 2021YFA1401900). Hongming Weng and Quansheng Wu were also supported by the Informatization Plan of the Chinese Academy of Sciences (Grant No. CASWX2021SF-0102). References Inhomogeneous Electron GasSelf-Consistent Equations Including Exchange and Correlation EffectsEfficient iterative schemes for ab initio total-energy calculations using a plane-wave basis setEfficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis setQUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materialsAdvanced capabilities for materials modelling with Quantum ESPRESSOFirst principles methods using CASTEPThe Abinitproject: Impact, environment and recent developmentsABINIT: Overview and focus on selected capabilitiesThe SIESTA method for ab initio order- N materials simulationS iesta : Recent developments and applicationsWIEN2k: An APW+lo program for calculating the properties of solidsDFTB+, a software package for efficient approximate density functional theory based atomistic simulationsEmpiric k·p Hamiltonian calculation of the band-to-band photon absorption in semiconductorsHidden Weyl points in centrosymmetric paramagnetic metalsMotion of Electrons and Holes in Perturbed Periodic FieldsImpact of Random Dopant Fluctuations on the Electronic Properties of In x Ga 1– x N/GaN Axial Nanowire HeterostructuresBand structure of indium antimonideTopological insulators from the perspective of first‐principles calculationsTopological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surfaceHexagonal Warping Effects in the Surface States of the Topological Insulator Bi 2 Te 3 Chern Semimetal and the Quantized Anomalous Hall Effect in HgCr 2 Se 4 Toward high-frequency operation of spin lasersThreshold current reduction in spin-polarized lasers: Role of strain and valence-band mixingSimulating the electronic properties of semiconductor nanostructures using multiband k · p models k · p theory for phosphorene: Effective g -factors, Landau levels, and excitonsCorrigendum: k.p theory for two-dimensional transition metal dichalcogenide semiconductors (2015 2D Mater. 2 022001)Ab Initio Studies of Exciton g Factors: Monolayer Transition Metal Dichalcogenides in Magnetic FieldsBand structure calculations of InP wurtzite/zinc-blende quantum wellsValley Zeeman effect and Landau levels in two-dimensional transition metal dichalcogenidesElectrons, holes, and excitons in GaAs polytype quantum dotsQuantum rings with Rashba spin-orbit coupling: A path-integral approachElectronic structure and Landé g-factor of a quantum ring in the presence of spin-orbit coupling: Temperature and Zeeman effectSpin–Orbit and Zeeman Effects on the Electronic Properties of Single Quantum Rings: Applied Magnetic Field and Topological DefectsHole energy levels and effective g-factor of quantum rings utilizing k.p Hamiltonian in terms of cylindrical polar coordinatesLandé g Factors and Orbital Momentum Quenching in Semiconductor Quantum DotsLande g-factor in semiconductor cylinder quantum dots under magnetic fields and spin-orbit interactionZeeman spin splittings in semiconductor nanostructuresElectron g factor in one- and zero-dimensional semiconductor nanostructuresOrbital Contributions to the Electron g Factor in Semiconductor NanowiresMesoscopic spin-orbit effect in the semiconductor nanostructure electron g factorLandé g Tensor in Semiconductor NanostructuresZeeman effect on surface electron transport in topological insulator Bi2 Se3 nanoribbonsTopological states and interplay between spin-orbit and Zeeman interactions in a spinful Su-Schrieffer-Heeger nanowireZeeman quantum-beat spectroscopy of NO2: Eigenstate-resolved Landé gF factors near dissociation thresholdPredicted Landé g-factors for open shell diatomic moleculesLandé g factors for 2p4 (3P)3p and 2p4(3P)3d states of Ne IIIrvsp: To obtain irreducible representations of electronic states in the VASPA k · p Effective Hamiltonian GeneratorMemorial Volume for Shoucheng ZhangLarge shift current, π Zak phase, and the unconventional nature of Se and TeIrRep: Symmetry eigenvalues and irreducible representations of ab initio band structuresDFT2kp: effective kp models from ab-initio dataProjector augmented-wave methodProjector augmented wave method:ab initio molecular dynamics with full wave functionsExponential protection of zero modes in Majorana islandsThe g ‐factor of the conduction electrons in Bi2 Se3Tunable effective g factor in InAs nanowire quantum dotsExciton g-factors in monolayer and bilayer WSe2 from experiment and theory
[1] Hohenberg P and Kohn W 1964 Phys. Rev. 136 B864
[2] Kohn W and Sham L J 1965 Phys. Rev. 140 A1133
[3] Kresse G and Furthmüller J 1996 Phys. Rev. B 54 11169
[4] Kresse G and Furthmüller J 1996 Comput. Mater. Sci. 6 15
[5] Giannozzi P, Baroni S, Bonini N, Calandra M, Car R, Cavazzoni C, Ceresoli D, Chiarotti G L, Cococcioni M, Dabo I, Corso A D, de Gironcoli S, Fabris S, Fratesi G, Gebauer R, Gerstmann U, Gougoussis C, Kokalj A, Lazzeri M, Martin-Samos L, Marzari N, Mauri F, Mazzarello R, Paolini S, Pasquarello A, Paulatto L, Sbraccia C, Scandolo S, Sclauzero G, Seitsonen A P, Smogunov A, Umari P, and Wentzcovitch R M 2009 J. Phys.: Condens. Matter 21 395502
[6] Giannozzi P, Andreussi O, Brumme T, Bunau O, Nardelli M B, Calandra M, Car R, Cavazzoni C, Ceresoli D, Cococcioni M, Colonna N, Carnimeo I, Corso A D, de Gironcoli S, Delugas P, DiStasio R A, Ferretti A, Floris A, Fratesi G, Fugallo G, Gebauer R, Gerstmann U, Giustino F, Gorni T, Jia J, Kawamura M, Ko H Y, Kokalj A, Küçükbenli E, Lazzeri M, Marsili M, Marzari N, Mauri F, Nguyen N L, Nguyen H V, Otero-de-la-Roza A, Paulatto L, Poncé S, Rocca D, Sabatini R, Santra B, Schlipf M, Seitsonen A P, Smogunov A, Timrov I, Thonhauser T, Umari P, Vast N, Wu X, and Baroni S 2017 J. Phys.: Condens. Matter 29 465901
[7] Clark S J, Segall M D, Pickard C J, Hasnip P J, Probert M I J, Refson K, and Payne M C 2005 Z. Kristallogr. 220 567
[8] Gonze X, Amadon B, Antonius G, Arnardi F, Baguet L, Beuken J M, Bieder J, Bottin F, Bouchet J, Bousquet E, Brouwer N, Bruneval F, Brunin G, Cavignac T, Charraud J B, Chen W, Côté M, Cottenier S, Denier J, Geneste G, Ghosez P, Giantomassi M, Gillet Y, Gingras O, Hamann D R, Hautier G, He X, Helbig N, Holzwarth N, Jia Y, Jollet F, Lafargue-Dit-Hauret W, Lejaeghere K, Marques M A L, Martin A, Martins C, Miranda H P C, Naccarato F, Persson K, Petretto G, Planes V, Pouillon Y, Prokhorenko S, Ricci F, Rignanese G M, Romero A H, Schmitt M M, Torrent M, van Setten M J, Troeye B V, Verstraete M J, Zérah G, and Zwanziger J W 2020 Comput. Phys. Commun. 248 107042
[9] Romero A H, Allan D C, Amadon B, Antonius G, Applencourt T, Baguet L, Bieder J, Bottin F C, Bouchet J, Bousquet E, Bruneval F, Brunin G, Caliste D, Côté M, Denier J, Dreyer C, Ghosez P, Giantomassi M, Gillet Y, Gingras O, Hamann D R, Hautier G, Jollet F C, Jomard G, Martin A, Miranda H P C, Naccarato F, Petretto G, Pike N A, Planes V, Prokhorenko S, Rangel T, Ricci F, Rignanese G M, Royo M, Stengel M, Torrent M, van Setten M J, Troeye B V, Verstraete M J, Wiktor J, Zwanziger J W, and Gonze X 2020 J. Chem. Phys. 152 124102
[10] Soler J M, Artacho E, Gale J D, García A, Junquera J, Ordejón P, and Sánchez-Portal D 2002 J. Phys.: Condens. Matter 14 2745
[11] García A, Papior N, Akhtar A, Artacho E, Blum V, Bosoni E, Brandimarte P, Brandbyge M, Cerdá J I, Corsetti F, Cuadrado R, Dikan V, Ferrer J, Gale J, García-Fernández P, García-Suárez V M, García S, Huhs G, Illera S, Korytár R, Koval P, Lebedeva I, Lin L, López-Tarifa P, Mayo S G, Mohr S, Ordejón P, Postnikov A, Pouillon Y, Pruneda M, Robles R, Sánchez-Portal D, Soler J M, Ullah R, Yu V W Z, and Junquera J 2020 J. Chem. Phys. 152 204108
[12] Blaha P, Schwarz K, Tran F, Laskowski R, Madsen G K H, and Marks L D 2020 J. Chem. Phys. 152 074101
[13] Hourahine B, Aradi B, Blum V, Bonafé F, Buccheri A, Camacho C, Cevallos C, Deshaye M Y, Dumitrică T, Dominguez A, Ehlert S, Elstner M, van der Heide T, Hermann J, Irle S, Kranz J J, Köhler C, Kowalczyk T, Kubař T, Lee I S, Lutsker V, Maurer R J, Min S K, Mitchell I, Negre C, Niehaus T A, Niklasson A M N, Page A J, Pecchia A, Penazzi G, Persson M P, Řezáč J, Sánchez C G, Sternberg M, Stöhr M, Stuckenberg F, Tkatchenko A, Yu V W Z, and Frauenheim T 2020 J. Chem. Phys. 152 124101
[14] Luque A, Panchak A, Mellor A, Vlasov A, Martí A, and Andreev V 2015 Physica B 456 82
[15] Gresch D, Wu Q, Winkler G W, and Soluyanov A A 2017 New J. Phys. 19 035001
[16] Luttinger J M and Kohn W 1955 Phys. Rev. 97 869
[17] Marquardt O, Geelhaar L, and Brandt O 2015 Nano Lett. 15 4289
[18] Kane E O 1957 J. Phys. Chem. Solids 1 249
[19] Zhang H J and Zhang S C 2013 Phys. Status Solidi RRL 7 72
[20] Zhang H J, Liu C X, Qi X L, Dai X, Fang Z, and Zhang S C 2009 Nat. Phys. 5 438
[21] Fu L 2009 Phys. Rev. Lett. 103 266801
[22] Xu G, Weng H M, Wang Z J, Dai X, and Fang Z 2011 Phys. Rev. Lett. 107 186806
[23] Faria J P E, Xu G F, Lee J, Gerhardt N C, Sipahi G M, and Žutić I 2015 Phys. Rev. B 92 075311
[24] Holub M and Jonker B T 2011 Phys. Rev. B 83 125309
[25] Marquardt O 2021 Comput. Mater. Sci. 194 110318
[26] Faria J P E, Kurpas M, Gmitra M, and Fabian J 2019 Phys. Rev. B 100 115203
[27] Kormányos A, Burkard G, Gmitra M, Fabian J, Zólyomi V, Drummond N D, and Fal'ko V 2015 2D Mater. 2 049501
[28] Deilmann T, Krüger P, and Rohlfing M 2020 Phys. Rev. Lett. 124 226402
[29] Faria J P and Sipahi G 2012 J. Appl. Phys. 112 103716
[30] Xuan F Y and Quek S Y 2020 Phys. Rev. Res. 2 033256
[31] Climente J I, Segarra C, Rajadell F, and Planelles J 2016 J. Appl. Phys. 119 125705
[32] Lucignano P, Giuliano D, and Tagliacozzo A 2007 Phys. Rev. B 76 045324
[33] Zamani A, Setareh F, Azargoshasb T, and Niknam E 2017 Superlattices Microstruct. 110 243
[34] León-González J C, Toscano-Negrette R G, Morales A L, Vinasco J A, Yücel M B, Sari H, Kasapoglu E, Sakiroglu S, Mora-Ramos M E, Restrepo R L, and Duque C A 2023 Nanomaterials 13 1461
[35] Zamani A and Rezaei G 2018 Superlattices Microstruct. 124 145
[36] Pryor C E and Flatté M E 2006 Phys. Rev. Lett. 96 026804
[37] Gharaati A 2017 Solid State Commun. 258 17
[38] Kotlyar R, Reinecke T L, Bayer M, and Forchel A 2001 Phys. Rev. B 63 085310
[39] Kiselev A A, Ivchenko E L, and Rössler U 1998 Phys. Rev. B 58 16353
[40] Winkler G W, Varjas D, Skolasinski R, Soluyanov A A, Troyer M, and Wimmer M 2017 Phys. Rev. Lett. 119 037701
[41] Toloza S M A, Ferreira D S A, de Andrada E S E A, and La R G C 2012 Phys. Rev. B 86 195302
[42] Alegre T P M, Hernández F G G, Pereira A L C, and Medeiros-Ribeiro G 2006 Phys. Rev. Lett. 97 236402
[43] Wang L X, Yan Y, Zhang L, Liao Z M, Wu H C, and Yu D P 2015 Nanoscale 7 16687
[44] Liu Z H, Entin-Wohlman O, Aharony A, You J Q, and Xu H Q 2021 Phys. Rev. B 104 085302
[45] Xin J and Reid S A 2002 J. Chem. Phys. 116 525
[46] Semenov M, Yurchenko S N, and Tennyson J 2016 J. Mol. Spectrosc. 330 57
[47] Fischer C F and Jönsson P 2001 J. Mol. Struct.: THEOCHEM 537 55
[48] Gao J C, Wu Q S, Persson C, and Wang Z J 2021 Comput. Phys. Commun. 261 107760
[49] Jiang Y, Fang Z, and Fang C 2021 Chin. Phys. Lett. 38 077104
[50] Song Z, Sun S, Xu Y, Nie S, Weng H, Fang Z, and Dai X 2021 First Principle Calculation of the Effective Zeeman's Couplings in Topological Materials (Singerpore: World Scientific) p 263
[51] Zhang R, Deng J, Sun Y, Fang Z, Guo Z, and Wang Z 2023 Phys. Rev. Res. 5 023142 The IR2PW code is available at https://github.com/zjwang11/IR2PW
[52] Iraola M, Mañes J L, Bradlyn B, Horton M K, Neupert T, Vergniory M G, and Tsirkin S S 2022 Comput. Phys. Commun. 272 108226
[53] Cassiano J A V V, Araújo A L, Junior P E F, and Ferreira G J 2023 arXiv:2306.08554 [cond-mat.mes-hall]
[54] Blöchl P E 1994 Phys. Rev. B 50 17953
[55] Blöchl P E, Först C J, and Schimpl J 2003 Bull. Mater. Sci. 26 33
[56] Albrecht S M, Higginbotham A P, Madsen M, Kuemmeth F, Jespersen T S, Nygård J, Krogstrup P, and Marcus C M 2016 Nature 531 206
[57]Zhang S, S H, Song Z D, Liang C, and Wang Z 2023 VASP2KP
[58] Köhler H and Wöchner E 1975 Phys. Status Solidi B 67 665
[59] Björk M T, Fuhrer A, Hansen A E, Larsson M W, Fröberg L E, and Samuelson L 2005 Phys. Rev. B 72 201307
[60] Förste J, Tepliakov N V, Kruchinin S Y, Lindlau J, Funk V, Förg M, Watanabe K, Taniguchi T, Baimuratov A S, and Högele A 2020 Nat. Commun. 11 4539