In the preparation step, we first identify the space group of the crystal structure for a nonmagnetic material with an even number of electrons per unit cell and standardize it using an open-source computational material science package, called Phonopy.^[30]

The structure data of a material may have numerical errors and different choices of coordinate axes and origin. We use Phonopy to standardize the structure to facilitate high-throughput calculations. Given an input structure and error tolerance, Phonopy will identify its corresponding space group, minimize the error of the atomic positions with respect to the space group, and transform the coordinate system to a consistent primitive cell basis convention. This consistent convention provides a basis to define the matrix forms of the symmetry operations in each space group that are used to calculate the irreps of the wavefunctions later.

We set a default error tolerance of 10⁻⁴ Å for Phonopy to identify the space group. If the input crystal structure is already associated with a space group and the designation is different from Phonopyʼs result, we increase the error tolerance until Phonopy finds the designated space group or the error tolerance reaches 0.1 Å, in which case the input data are considered to have large numerical errors and should be replaced by a more accurate structure.

Take SnTe of space group 225 as an example. Table 1 lists its primitive cell basis vectors of the original and standardized structures. Notice that the coordinate system is transformed.

Table 1.

Table 1.

Table 1. Primitive cell basis vectors of SnTe in the soc setting (unit: Å). . Original structure (3.8690, 0.0000, 2.2338) (1.2897, 3.6477, 2.2338) (0.0000, 0.0000, 4.4675) Standardized structure (0.0000, 3.1590, 3.1590) (3.1590, 0.0000, 3.1590) (3.1590, 3.1590, 0.0000)

Table 1. Primitive cell basis vectors of SnTe in the soc setting (unit: Å). .

The standardized crystal structure is fed to the first-principles calculation software once obtained. Following the standard DFT calculation procedure, the self-consistent electron density is first calculated and used in the second step to compute the wavefunctions along with the energy eigenvalues at a given list of HSPs in the BZ. The lowest N bands (N being the number of electrons per unit cell) at each HSP are defined as the valence bands.

An HSP is a k-point in the BZ with the highest symmetry in its neighborhood. The positions of the HSPs are determined by the space group with a maximum number of 8 in all 230 space groups. The coordinates of the HSPs for each space group can be found at “The k-vector types and Brillouin zones of Space Groups” on the Bilbao Crystallographic Server (BCS).^[31,32] Two HSPs related by time-reversal symmetry (with opposite coordinates and one with a suffix “A,” e.g., “K” and “KA”) are both listed on this page, and we only use one of them.

The HSPs and the symmetry operations of the 230 space groups on the BCS are given in its conventional cell basis convention. They must be transformed into Phonopyʼs primitive cell basis convention to be consistent with the crystal structure and the DFT calculation results. These two sets of convention differ by a similarity transformation of the coordinate axes and a shift of origin, which are listed in Appendix B. In either case, the coordinates of the HSPs are given as fractions of the reciprocal lattice basis, with the reciprocal lattice basis vector being defined by the standard formulas

SnTe of space group 225 only has four HSPs in the BZ, namely , L = [1/2, 1/2, 1/2], W = [1/2, 1/4, 3/4], and X = [1/2, 0, 1/2], with the fractional number being the coordinates under the reciprocal lattice basis in Phonopyʼs primitive cell convention.

We conducted the calculation with the generalized gradient approximation of the Perdew–Burke–Ernzerhof (PBE)-type exchange-correlation potential using the pseudopotential files of the projector augmented-wave atomic data (PAW_GGA_PBE). Subsection 3.1 lists the other important input parameters. The calculation may be performed with the spin–orbit coupling turned on and off, which are briefly called the “soc setting” and the “nsoc setting,” respectively. Note that although spin–orbit coupling always exists in realistic materials, ignoring it is a more physically relevant approximation for systems consisting of small atoms. For materials with a large spin–orbit coupling, the results obtained in the nsoc setting can be used as a reference and may help diagnose the origin of topological properties in the soc setting.

Extracting the irrep of each (multiplet of) valence band(s) is one of the core modules of the algorithm. In the case of one-dimensional irreps (the higher dimensional case will be discussed later), the wavefunction of each band at each HSP has the corresponding little group symmetry and is a basis function of one of the little groupʼs irreps. For such a wavefunction, its corresponding irrep is extracted in two steps. First, we calculate its character with each symmetry operation in the little group. Second, we compare these characters to the character tables on the BCS to identify the irreps. Our algorithm for extracting irreps currently only supports wavefunctions in a plane wave basis. We discuss the main steps in the following paragraphs.

Assume that the wavefunction in the nsoc setting is expanded in the plane wave basis:

In the case of one-dimensional irreps, the character of under a certain symmetry operation is simply the eigenvalue λ :

However, the representation of the translation subgroup of each space group on BCS takes the convention

We identify the irrep of the wavefunction when the characters of all the symmetry operations are calculated. This step uses the orthogonal theorem in group representation theory:

Note that although a set of bands belonging to a single high-dimensional irrep is, in principle, degenerate, a small energy difference between such bands may be observed because of numerical errors. As a result, we set a degeneracy error of min (0.5 meV, , where is the direct band gap at the k-point. We split the valence bands into small sets, in which bands have an energy difference smaller than this error, and calculate the irrep for each set of bands separately.

Moreover, time-reversal symmetry always exists in nonmagnetic materials, thereby causing an additional degeneracy between the original irreps of the space groups. An irrep at a high-symmetry momentum is generally mapped under time reversal to either (i) itself or (ii) another different irrep or (iii) another same irrep. As a result, time-reversal symmetry puts different constraints on the occupation numbers of different types of irreps in the valence bands: For type-(i) irreps, the occupation numbers do not have additional constraints; for type-(ii) irreps, the occupation numbers for the two irreps related by time-reversal symmetry are constrained to equal; and for type-(iii) irreps, the occupation number of a type-(iii) irrep must be even. The compatibility relations are given in terms of the linearly independent occupation numbers; hence, we must extract them after considering these constraints.

The symmetry operations and the character tables of the 230 space groups are obtained from the “Irreducible representations of the double space groups” database on the BCS. Note that two tables can be found on the page, and we only use the first one (i.e., representation matrices of the little group). The additional degeneracy caused by the time-reversal symmetry is not included in this page, and we extract this information separately from the “Band representations and elementary band representations of double space groups” database on the BCS, in which the irreps written together like “ ” are degeneracies caused by time-reversal symmetry (the number 4 in the parentheses denotes the dimension of the whole degenerate irrep).

As mentioned in Subsection 2.2, the symmetry operations on the BCS are given in its conventional cell convention. We transform their coordinate axes and origin to make them consistent with Phonopyʼs primitive cell convention. A change of the coordinate axes transforms the O(3) matrix of the symmetry operation, whereas a different choice of origin changes the translation part by the formula

The SU(2) part of the symmetry operations is directly calculated using the standard formula as follows:

Table 2 lists the calculated irreps of SnTe in the soc setting, with the number in the parentheses denoting the irrep dimension. Notice that SnTe has 20 electrons per unit cell (i.e., it has 20 valence bands); this number of electrons equals the sum of the number of irreps weighted based on their dimensions at each HSP. The time-reversal symmetry will pin and together. SnTe has the same number of these two irreps.

Table 2.

Table 2.

Table 2. Irreps of SnTe in the soc setting. . Irreps of the (2) (2) (2) (4) (4) number 2 1 1 2 1 Irreps of the L (1) (1) (2) (2) number 3 3 5 2 Irreps of the W (2) (2) number 5 5 Irreps of the X (2) (2) (2) (2) number 4 3 2 1

Table 2. Irreps of SnTe in the soc setting. .

The procedure for extracting the band irreps described earlier is applicable to both valence and conduction bands. In our algorithm, we first calculate the irreps of the bands near the top of the valence band after obtaining the wavefunctions and energy eigenvalues at the HSPs from the DFT calculations to determine if the material is an HSPSM.

A material with partially filled irreps at some HSPs is labeled as an HSPSM. In other words, the material has degeneracy between valence and conduction bands at some HSPs, and this degeneracy is protected by the HSPʼs symmetry, making the material a robust topological semimetal. Note that the degeneracy can be caused by space group symmetry and time-reversal symmetry.

In principle, an HSPSM is gapless at each degenerate HSP, but a small direct gap may exist because of the numerical errors of the first-principle calculations. We perform an HSPSM examination for all materials with a direct gap less than 2 meV. An energy tolerance for degeneracy starting from 0 is set starting from the N-th band (N being the number of electrons per unit cell). The lower and upper bands of the N-th band within this tolerance are found, and their irreps are calculated. If the integer number of the irreps is successfully identified, we check whether the N-th and (N + 1)-th band belong to the same irrep. If so, the material is labeled as an HSPSM; otherwise, we proceed to the next phase of the algorithm. By contrast, if we fail to obtain the integer number of the irreps, the energy tolerance is increased to include more bands. If this tolerance reaches 2 meV, and the integer number of the irreps still cannot be identified, the material is considered to be poorly converged and should be re-calculated in the DFT phase using a higher-convergence threshold. The highest dimension of irreps in all the space groups is 8; therefore, we only need to examine at most eight bands in this step.

Note that in principle, only the information on the valence bands is required to determine the HSPSMs because the of the top valence band is always degenerate with at least one conduction band in an HSPSM if (i) the linear decomposition of the valence bands into the irreps has a non-integer coefficient or if (ii) the coefficients of the irreps that are related to each other (i.e., the type-(ii) irreps) are not equal or if (iii) the coefficient of a type-(iii) irrep is odd. These factors can be used as the criteria for determining the HSPSMs (i.e., whether all the valence bands form an integer number of symmetry-enforced degenerate multiplets). However, we also include several low-conduction bands in our algorithm to obtain the degeneracy dimension.

The direct gaps for SnTe at the four HSPs in the soc setting are 4.4196, 0.0939, 1.6635, and 5.6673 eV, with the smallest one being larger than 2 meV, thereby excluding the possibility for SnTe to be an HSPSM.

Take HgTe of space group 216 as another example. It has a zero direct gap at in the soc setting and has 18 electrons per unit cell. Our algorithm confirms that the 17th, 18th, 19th, and 20th bands at belong to the same four-dimensional irrep , making HgTe an HSPSM in the soc setting.

We continue to extract the irreps of all the valence bands and check them against the compatibility conditions on each HSL if a material has no degeneracy between the valence and conduction bands at any HSP (i.e., it is not an HSPSM). The material is considered to be a band insulator and the algorithm proceeds into the next step if all the compatibility conditions are satisfied; otherwise, the material is labeled as an HSLSM, indicating that symmetry-enforced band crossings exist between the conduction and valence bands along the HSLs, where the compatibility conditions are violated.

The little group of an HSL is a subgroup of the little group of any HSP on it; this means that each irrep of the little group at an HSP can be decomposed into the direct sum of the irreps of the little group at a joining HSL. The decomposition coefficients are described by the so-called compatibility relations. For a system to be gapped along an HSL AB that connects two HSPs A and B, the valence bands of A and B must be continuous; thus, their irreps at A and B must be decomposable into the same set of irreps of ABʼs little group with exactly the same coefficients. We call such a condition the compatibility condition, which can be represented as a set of linear equations involving the total number of each irrep for all valence bands at A and B. The compatibility condition being satisfied along all the HSLs is a necessary condition for a system to have a direct gap. Note that if a material satisfies all the compatibility conditions, band crossings may still exist at some HSLs between the valence and conduction bands; however, the crossing point(s) are not symmetry-enforced and may be gapped by tuning the Hamiltonian without breaking the symmetry. Moreover, note that the HSL mentioned herein is not necessarily a line. It can also be a surface. The generic points in the BZ are also named as an HSL “GP” for convenience, which only gives trivial compatibility conditions (i.e., conservation of the number of bands between two k points).

The compatibility relations for the 230 space groups have been explicitly derived and are available at the “Compatibility Relations between representations of the Double Space Groups” database on the BCS. For each space group, the BCS lists these relations between each HSP/HSL and the HSP/HSL connected to it. We use these data to construct two matrices of compatibility conditions for each HSL, that is, for the soc and nsoc settings, respectively. Element (i, j) of the matrix gives the decomposition coefficient of the j-th HSP irrep into the i-th HSL irrep, with a minus sign if the irrep belongs to the second HSP. The irrep vector of a material is generated by collecting the irreps of all valence bands at both HSPs, with each entry being the total number of occurrences of one irrep arranged in the same order matching the columns of the compatibility condition matrix. We obtain a vector, whose nonzero elements indicate the violated compatibility conditions, by multiplying the compatibility condition matrix to the irrep vector.

As an example, Figure 2 illustrates one page of the compatibility relations of space group 225 on the BCS, which lists the compatibility relations between HSL , ) and the three HSP/HSLs connected to it, namely , W = (1/2, 1, 0), and generic points GP.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Compatibility relations involving the k-vector Q (No. 225) on the BCS.

The irreps , , , and W₁, W₄, W₅ for L and W are connected to Q₁, whereas , , , and W₂, W₃, W₅ are connected to Q₂. These two relations can be transformed into two linear equations, namely

They can be written in a matrix form, as presented in Table 3.

Table 3.

Table 3.

Table 3. Compatibility condition matrix between L and W. . Irreps (1) (1) (1) (1) (2) (2) W1(1) W2(1) W3(1) W4(1) W5(2) 1 1 0 0 1 1 −1 0 0 −1 −1 0 0 1 1 1 1 0 −1 −1 0 −1

Table 3. Compatibility condition matrix between L and W. .

The whole compatibility condition matrix can be constructed similarly, with each row being one such linear equation where coefficients for all irrelevant irreps are set to 0.

The irreps of SnTe satisfy all the compatibility conditions in the soc setting but not in the nsoc setting. Table 4 lists the irreps of SnTe at L and W in the nsoc setting.

Table 4.

Table 4.

Table 4. Irreps of SnTe in the nsoc setting. . Irreps (1) (1) (1) (1) (2) (2) W1(1) W2(1) W3(1) W4(1) W5(2) Number 2 0 0 2 3 0 3 2 0 1 2

Table 4. Irreps of SnTe in the nsoc setting. .

They violate both the L–Q₁–W and L–Q₂–W compatibility conditions, making SnTe an HSLSM in the nsoc setting.

We calculate the symmetry-based indicators for a material if its band structure satisfies all the compatibility conditions of the space group. The nonzero calculated indicators are labeled as a GMSM in the nsoc setting or a TI or TCI in the soc setting. The parity of the last number in the indicators in the soc setting determines whether the material is a TI or a TCI, with an odd number corresponding to a TI and an even one to a TCI. Note that not all space groups have symmetry-based indicators (117 space groups in the soc setting and 53 in the nsoc setting). A material is labeled as a trivial insulator if it has zero-valued indicators or belongs to a space group with no indicator.

The group structures of the symmetry-based indicators are given in Ref. [14], and their explicit formulas in terms of the corresponding space group irreps used in our algorithm are derived in Refs. [15] and [18]. As demonstrated by Ref. [18], the nonzero indicators in the nsoc setting correspond to the topological semimetals with topological nodes at the non-high-symmetry momenta. The possible configurations of these nodes can be found in the same article. The mapping between the symmetry indicators and the possible sets of topological invariants is tabulated in Ref. [15]. The mapping is generally one-to-many, with one set of indicators mapping to several possible nonequivalent sets of topological invariants.

Note that materials with zero indicators or belonging to the space groups with no indicators are not necessarily topologically trivial. This only indicates that their band topology, if existing, cannot be detected using any symmetry eigenvalue diagnosis approach because considering only the symmetry information at the HSPs is naturally incomplete.

The symmetry-based indicator formula of space group 225 in the soc setting is^[15]

The SnTe indicator is calculated as , indicating that SnTe is a TCI, by substituting the number of irreps of the valence bands of SnTe listed in Subsection 2.3.

The input generation module standardizes the structure file provided by the user, as mentioned in Subsection 2.1, and prepares all the necessary input files for VASP or ABINIT to run the DFT calculation. Several computational material science packages, including Pymatgen,^[33] AbiPy^[28,34] (only for ABINIT), and Phonopy, are used by the module to accomplish the task.

The DFT calculation involves a self-consistent electron density calculation and a subsequent step to compute the wavefunctions and eigenvalues at all relevant HSPs. For VASP, the input files for these two steps are stored separately in two folders named “scf” and “wave,” with each folder containing four files “POSCAR,” “INCAR,” “KPOINTS,” and “POTCAR,” which store the structure, parameter, k points, and pseudopotential information in a VASP-native format, respectively.

For ABINIT, both steps are defined in the same set of input files, including a “.in” file that summarizes the structure, parameter, and k-point information and a “.files” file that designates the name and the location of the pseudopotential files, the “.in” file, and several output files.

Table 5 presents a summary of the important default input parameters for VASP and ABINIT assumed by the input generation module. These default setup yields satisfying results in most cases. Users can also manually adjust these input files to tune the parameters for the DFT calculation when necessary. For example, in case the first-principle calculations do not converge, users could retry the calculation with a different integral method for the k points or a different convergence scheme.

Table 5.

Table 5.

Table 5. Key input parameter template for VASP and ABINIT. . Software VASP ABINIT Common parameters Cut-off energy PREC = higha ecut and pawecutdg are set according to the pseudopotential file Convergence criterion NELM = 300 nstep = 300 EDIFF = 1e-06 toldfe = 1.0d — 6 eV tolwfr2 = 1.0d — 10 eV Soc-setting LSORBIT = True pawspnorb = 1 MAGMOM = 3*N_atom*0.0 magconon = 1 Nsoc-setting LSORBIT = False pawspnorb = 0 MAGMOM = N_atom*0.0 magconon = 1 istwfk2 = num_k*1 Other parameters ISTART = 0 iomode = 3 LMAXMIX = 4 fband = 1 Parameters for self-consistent calculation Occupation option ISMEAR = −5 occopt1 = 3 Charge density ICHARG = 2 getden1 = 0 prtden1 = 1 Wave function LWAVE = False prtwf1 = −1 K-points Auto nshiftk1 = 1 30 shiftk1 = 0.0 0.0 0.0 (in KPOINTS file) ngkpt1 = 3 3 3 kptopt1 = 1(nsoc)/4(soc) Parameters for band calculation Occupation option ISMEAR = 0 occopt2 = 3 Charge density ICHARG = 11 getden2 = −1 Wave function LWAVE = TRUE prtwf2 = 1 prteig2 = 1 enunit2 = 1 K-points List the coordinates of HSPs (in KPOINTS file) nkpt2 = number of k points kpt2 = coordinates of HSPs a This setting automatically sets ENCUT = 125% ENMAX, where ENMAX is read from the pseudopotential file.

Table 5. Key input parameter template for VASP and ABINIT. .

The module currently works with the PAW_GGA_PBE-type pseudopotential files. Users need to configure the package with the location of the proper pseudopotential files before using the module. In using the other types of pseudopotential files, the user could use the structure and k-point information automatically generated by the module and manually edit the input files to fill in the correct values for the other parameters.

The analysis module of the package takes in the output files of the first-principle calculation and reports the topological classification of the material. For VASP, the output files from the band calculation, including “OUTCAR,” “WAVECAR,” and “EIGENVAL,” are analyzed. For ABINIT, the “_WFC.nc” file that stores the wavefunction in NetCDF format is analyzed.

Aside from a classification, the module also reports relevant details for the topological class. For an HSPSM, the module lists the HSPs and the corresponding irreps, where degeneracy exists between the valence and conduction bands. For an HSLSM, the module lists the pairs of HSPs, where the band crossings between the valence and conduction bands exist at some HSLs connecting these two HSPs, and the irreps of all valence bands at each HSP. For a GMSM/TI/TCI, the symmetry-based indicators and the irreps of all valence bands are given. Lastly, only the triviality and band irreps are listed for a trivial insulator.