qeirreps: an open-source program for Quantum ESPRESSO to compute irreducible representations of Bloch wavefunctions

Let us consider a space group $\mathcal{G}$. An element $g\in\mathcal{G}$ moves $\bm{x}\in\mathbb{R}^{3}$ to

where $p_{g}\in\text{O}(3)$ represents the point group part and $\bm{t}_{g}\in\mathbb{R}^{3}$ represents the translation part of $g$. The product of two elements $g\in\mathcal{G}$ and $g^{\prime}\in\mathcal{G}$ is defined by

In general, $\bm{t}_{g}$ is not necessarily a lattice vector. If we can choose the origin of $\bm{x}$ in such a way that $\bm{t}_{g}$ becomes a lattice vector simultaneously for all $\forall g\in\mathcal{G}$, the space group is symmorphic; otherwise, it is nonsymmorphic. Nonsymmorphic space groups typically contain either screw axes or glide planes.

A space group $\mathcal{G}$ always has a subgroup $T$ composed of lattice translations. An element $T_{\bm{R}}$ of $T$ can be characterized by a lattice vector $\bm{R}$. We introduce the wavevector $\bm{k}$ through the representation of the translation symmetry $U_{\bm{k}}(T_{\bm{R}})=e^{-i\bm{k}\cdot\bm{R}}$.

An element $g\in\mathcal{G}$ changes $\bm{k}$ to $g(\bm{k})=p_{g}\bm{k}$. We say $\bm{k}$ is invariant under $g$ when $g(\bm{k})=\bm{k}+\bm{G}$ with a reciprocal lattice vector $\bm{G}$. For each $\bm{k}$, the set of $g\in\mathcal{G}$ that leaves $\bm{k}$ invariant, i.e., $\left\{h\in\mathcal{G}\ |\ h(\bm{k})=\bm{k}+\bm{G}\right\}$, forms a subgroup of $\mathcal{G}$ called the little group $\mathcal{G}_{\bm{k}}$. Since $\mathcal{G}_{\bm{k}}$ is also an infinite group due to its translation subgroup $T$, sometimes the finite group $\mathcal{G}_{\bm{k}}/T$, called “little co-group [72],” is discussed instead. In this work, we always consider the little group $\mathcal{G}_{\bm{k}}$.

The space group symmetry $\mathcal{G}$ can be encoded in the single-particle Hamiltonian $H_{\bm{k}}$ in momentum space by requiring

for each $g\in\mathcal{G}$, where $U_{\bm{k}}(g)$ is a unitary matrix forming a representation of $\mathcal{G}$. For spinful electrons, relevant representations become ‘projective’ (in contrast to the ordinary ‘linear’ representation) and $U_{\bm{k}}(g)$’s satisfy

for $g,g^{\prime}\in\mathcal{G}$. The projective factor $\omega^{\text{sp}}(g,g^{\prime})=\pm 1$ is not unique, and we choose it as in Eq. (9). For spin-rotation invariant systems, $\omega^{\text{sp}}(g,g^{\prime})$ can be set $1$ by neglecting the spin degree of freedom.

There is a way of treating projective representations of $\mathcal{G}$ as linear representations of an enlarged group $\mathcal{G}^{\prime}$ who has the doubled number of elements. This method is known as ‘double group,’ and is commonly used, for example, in Ref. [73]. In this treatment, $\mathcal{G}^{\prime}$ contains both $\pm g$ for each element of $\mathcal{G}$, and the product of $\eta g\in\mathcal{G}^{\prime}$ and $\eta^{\prime}g^{\prime}\in\mathcal{G}^{\prime}$ is defined by $\eta\eta^{\prime}\omega^{\text{sp}}(g,g^{\prime})gg^{\prime}\in\mathcal{G}^{\prime}$. See C for an example. Whether one handles projective representations directly, as we do in our program, or via the double group method does not affect the physical conclusion. See more detailed discussions on the relationship between projective representations and double groups in Refs. [74, 75].

The single-particle Hamiltonian $H_{\bm{k}}$ can be diagonalized as

where $\epsilon_{\bm{k}n}$ ($n=1,2,\cdots,M_{\bm{k}}$) is the $n$-th energy level of $H_{\bm{k}}$ and $\mathds{1}_{n}$ is the identity matrix whose size is given by the order of degeneracy of $\epsilon_{\bm{k}n}$. The unitary matrix $\Psi_{\bm{k}}$ is composed of all eigenvectors of $H_{\bm{k}}$. When $h$ is an element of $\mathcal{G}_{\bm{k}}$, $U_{\bm{k}}(h)$ can also be block-diagonalized by $\Psi_{\bm{k}}$ as

Here, $U_{\bm{k}n}(h)$ is the representation of $\mathcal{G}_{\bm{k}}$ of the $n$-th band. Although the specific representation depends on the detailed choice of $\Psi_{\bm{k}}$, its character

is basis independent. The output of our program is the list of characters $\chi_{\bm{k}n}(h)$ for each the energy level $\epsilon_{\bm{k}n}$ at each high-symmetry point $\bm{k}$.

In the absence of accidental degeneracy or symmetries other than $\mathcal{G}$, the representation $U_{\bm{k}n}(h)$ is automatically irreducible. That is, the character $\chi_{\bm{k}n}(h)$ of the $n$-th band coincides with one of the characters $\chi_{\bm{k}}^{\alpha}(h)=\mathrm{tr}\left[U_{\bm{k}}^{\alpha}(h)\right]$ of irreducible representations of $\mathcal{G}_{\bm{k}}$. Otherwise, $U_{\bm{k}n}(h)$ can be decomposed into irreducible representations $U_{\bm{k}}^{\alpha}$ as $U_{\bm{k}n}=\oplus_{\alpha}n_{\bm{k}n}^{\alpha}U_{\bm{k}}^{\alpha}$, where the multiplicity of each irreducible representation is given by the formula

Our output file provides the list of $\chi_{\bm{k}}^{\alpha}(h)$ and $n_{\bm{k}n}^{\alpha}$ in irreps_list.dat and irreps_number.dat, respectively. ¹¹1Note that the output of our program is the character of representations of $\mathcal{G}_{\bm{k}}$, not $\mathcal{G}_{\bm{k}}/T$, and one should not confuse them. Irreducible representations of $\mathcal{G}_{\bm{k}}$ and $\mathcal{G}_{\bm{k}}/T$ are related to each other by a simple rule, and one can convert one to the other easily.

The sum of integers $n_{\bm{k}n}^{\alpha}$ in Eq. (8) over all filled bands, i.e., $n_{\bm{k}}^{\alpha}=\sum_{n:\text{filled}}n_{\bm{k}n}^{\alpha}$, counts the irreducible representation $U_{\bm{k}}^{\alpha}(h)$ below the Fermi level at each high-symmetry point $\bm{k}$. The integers $n_{\bm{k}}^{\alpha}$ can be used, for example, to diagnose the topology of band insulators through the method of symmetry indicators [51, 52] or the topological quantum chemistry [53]. These methods diagnose the topology of the target material by comparing its irreducible representations with those of atomic insulators. Useful formulas of topological indices in terms of $\{n_{\bm{k}}^{\alpha}\}$ are provided in Refs. [76, 77, 78, 60].

For the user’s convenience, we implemented the function that automatically computes the $\mathbb{Z}_{4}$ index as the sum of the inversion parities of all occupied bands at all the time-reversal invariant momenta (TRIMs). The $\mathbb{Z}_{4}$ index can detect not only strong topological insulators but also topological crystalline insulators such as mirror Chern insulators or higher-order topological insulators [77, 78]. See Sec. 4.1 for an example. Other symmetry indicators can also be computed in the same way using the output of our program. We discuss several examples in Sec. 4.4.

Let us summarize our conventions that are necessary to compare the output of our program with that of others. There are three sources of ambiguities that affect the U(1) phase of the representation $U_{\bm{k}}^{\alpha}(h)$: (i) the choice of representatives of $\mathcal{G}$, (ii) the choice of the coordinates, and (iii) the choice of spin rotation matrices. Readers not interested in the details can skip to Sec. 3.

A Fortran 90 compiler, BLAS, and LAPACK libraries are required for the installation of qeirreps. Quantum ESPRESSO must also be installed in advance. In addition, the program qe2respack [79] ²²2This program qe2respack is a part of the program package RESPACK [79]. qe2respack in the latest version of RESPACK does not support DFT calculations with spin-orbit coupling. For calculations with spin-orbit coupling, we need the specific version of qe2respack described in this section. is required, which can be downloaded or cloned from the branch of respack “maxime2” in the GitHub repository (https://github.com/mnmpdadish/respackDev/). The program qe2respack.py is in the directory util/qe2respack.

Our program qeirreps is released at GitHub (https://github.com/qeirreps/qeirreps). The program files which contain source files, documents, and examples can be cloned or downloaded from this repository. The file qeirreps/src/Makefile must be edited to specify the compiler and libraries in your compiling environment. By typing $ make in the source directory qeirreps/src/, the executable binary qeirreps.x is compiled.

Our program qeirreps works based on the output of Quantum ESPRESSO. To prepare input files of qeirreps, the following three steps are needed to be done one by one:

1. 1.

   Self-consistent first-principles (scf) calculation of a target material.

2. 2.

   Non-self-consistent first-principles (nscf) calculation of the material for each high-symmetry momentum.

3. 3.

   Data conversion from Quantum ESPRESSO output files to qeirreps input files.

The first two steps can be carried out by the standard functions of Quantum ESPRESSO. The result of the scf calculation is used in the successive nscf calculation. The set of high-symmetry points to be included in the calculation depends on the purpose of the calculation and the space group of the target material. When one applies the results of qeirreps to symmetry indicators or topological quantum chemistry, the minimum set of high-symmetry points are listed in Refs. [78, 57, 58].

We refer to the directory that stores the wavefunction data computed by the nscf calculation OUTDIR/PREFIX.save below. For the moment, qeirreps requires norm-conserving calculations. Pseudo-potentials must be optimized for norm-conserving calculations (for example, Optimized Norm-Conserving Vanderbilt Pseudopotential [80] is available in PseudoDojo (http://www.pseudo-dojo.org) [81]) and the option wf_collect should be set .TRUE. in the input files for DFT calculations.

The output files in OUTDIR/PREFIX.save must be converted by qe2respack into the form of input files of qeirreps. The work directory is referred to as DIRECTORY_NAME here. Then, type the following in the work directory:

where PATH_OF_qe2respack is the path to the directory util/qe2respack. qeirreps reads the files produced by qe2respack.py in the dir-wfn directory. The contents of these files are summarized in Table 1.

To run qeirreps, a directory named output must be created in the directory DIRECTORY_NAME above. Then type

where PATH_OF_qeirreps is the path to the directory qeirreps/src. There should be 12 output files in output. Seven of them named *.dat contain the following information: (i) the list of the point group part $p_{g}$ in pg.dat, (ii) the list of the translation part $\bm{t}_{g}$ in tg.dat, (iii) the list of the spin rotation part $p^{sp}_{g}$ in srg.dat, (iv) the factor system $\omega^{\text{sp}}(g,g^{\prime})$ associated with the electronic spin in factor_system_spin.dat, (v) the representation characters $\chi_{\bm{k}n}(h)$ of Bloch wavefunctions in character.dat, (vi) the characters $\chi_{\bm{k}}^{\alpha}$ of irreducible representations of $\mathcal{G}_{\bm{k}}$ in irreps_list.dat, and (vii) the numbers $n_{\bm{k}n}^{\alpha}$ of irreducible representations of $\mathcal{G}_{\bm{k}}$ in the $n$-th energy level in irreps_number.dat. The contents of these files are summarized in Table 2. Five of them named *_import.txt are the corresponding files for Mathematica usage. The standard output that appears during this calculation shows the lattice vectors, reciprocal lattice vectors, operation type of each symmetry operators, and so on.

For materials with inversion symmetry, qeirreps also implements the automatic evaluation of the $\mathbb{Z}_{4}$ index. To do this, an option of filling should be added to the command as

Here, z4 is the option to obtain the $\mathbb{Z}_{4}$ index and FILLING is the number of electrons per unit cell of the target material and is shown in the standard output of scf calculation by Quantum ESPRESSO as “number of electrons = FILLING.” Then, the program generates an additional output named z4.dat that contains the value of the $\mathbb{Z}_{4}$ index.

Our program qeirreps also provides the input file for the tool CheckTopologicalMat (https://www.cryst.ehu.es/cgi-bin/cryst/programs/topological.pl) [57]. Our program exports the input file named Bilbao.txt with the following command

where ctm is the option to obtain the input file. The tool CheckTopologicalMat tells us topological properties such as the value of symmetry indicators and gaplessness.

Our first example is bismuth. The space group ${R\bar{3}m}$ (No.$166$) contains the inversion symmetry $I$. We obtain the crystal structure data of bismuth [Figure 3 (a)] from Material Project [82] and converted it into the form of an input file for Quantum ESPRESSO (included in A.1) by SeeK-path [83, 84] ³³3In this particular example, we manually shifted positions of atoms in the primitive unit cell so that we can compare our results with the one in Ref. [57].

We compute the irreducible representations of wavefunctions by qeirreps, taking into account the spin-orbit coupling. Our results are shown in the band structure in Figure 3 (b). The correspondence between the labels and characters of irreducible representations is included in B.1. These results are consistent with previous studies [57, 73].

We also compute the $\mathbb{Z}_{4}$ index [77, 78] using the option explained in Sec. 3.4. The output (z4.dat) shows

Namely, the $\mathbb{Z}_{4}$ index for bismuth with significant spin-orbit coupling is $2$, indicating that this material is a higher-order topological insulator [85, 58, 57, 86]. The sample input files for bismuth are available in the directory qeirreps/example.

To demonstrate that qeirreps equally works for nonsymmorphic space groups, let us discuss silicon. Its space group is ${Fd\bar{3}m}$ (No. $227$), which also contains the inversion symmetry $I$. The calculation procedure is completely the same as those for bismuth in Sec. 4.1. Our results of irreducible representations are in Figure 4 (b). The $\mathbb{Z}_{4}$ index is found to be zero. These results are consistent with the previous study in Ref. [57]. Sample input files for silicon are also in the directory qeirreps/example.

Here we discuss NaCdAs as an example of topological materials with nonsymmorphic space group symmetries. Its space group is $Pnma$ (No. $62$) which has the inversion symmetry $I$. Irreducible representations computed in the presence of spin-orbit coupling are shown in Figure 5 (b). The $\mathbb{Z}_{4}$ index is 1, indicating that this material is a candidate of the strong $\mathbb{Z}_{2}$ topological insulator. These results are consistent with the previous study in Ref. [57]. (The $\mathbb{Z}_{4}$ index in Ref. [57] is 3 because their definition of the index includes an additional minus sign.) Sample input files for NaCdAs are contained in the directory qeirreps/example.

As our fourth example, let us discuss PbPt₃, whose space group ${Pm\bar{3}m}$ (No. $221$). This space group contains various elements such as the inversion $(I)$, the rotoinversion about the $z$-axis $(S_{4}^{z})$, and the mirror symmetries. We also assume the time-reversal symmetry. Our results of irreducible representations in the presence of spin-orbit coupling are in Figure 6 (b).

In this symmetry settings, we can define a $\mathbb{Z}_{4}$ index (other than the sum of the inversion parities) and a $\mathbb{Z}_{8}$ index by [87, 77, 78]

where $n_{K}^{\alpha_{K}}$ represents the number of occupied states that have irreducible representations $U_{K}^{\alpha_{K}}$ of $\mathcal{G}_{K}$, and the set of $\{n_{K}^{\alpha_{K}}\}$ is included in B.8. $K_{4}$ denotes four momenta invariant under the $S_{4}^{z}$. Although our program qeirreps does not offer an automated calculation of these indices, one can compute them manually by using character.dat or uploading Bilbao.txt to CheckTopologicalMat [57]. In this example, we find $(z_{4},z_{8})=(3,6)$, which is consistent with Ref. [57]. According to Ref. [77, 78], these values indicate nontrivial mirror Chern numbers. Sample input files for PbPt₃ are included in qeirreps/example.