WannierTools: An open-source software package for novel topological materials

TL;DR

WannierTools is an open-source software for analyzing topological materials, capable of classifying phases, identifying Weyl/Dirac points, and calculating Berry phases and surface states within the tight-binding framework.

Contribution

It introduces a comprehensive tool that integrates various topological analysis functions for materials, enhancing research capabilities in topological condensed matter physics.

Findings

Classifies topological phases using Wilson loop calculations

Identifies Weyl/Dirac points and nodal lines in materials

Calculates Berry curvature and surface state spectra

Abstract

We present an open-source software package WannierTools, a tool for investigation of novel topological materials. This code works in the tight-binding framework, which can be generated by another software package Wannier90. It can help to classify the topological phase of given materials by calculating the Wilson loop and can get the surface state spectrum which is detected by angle-resolved photoemission (ARPES) and in scanning tunneling microscopy (STM) experiments. It also identifies positions of Weyl/Dirac points and nodal line structures, calculates the Berry phase around a closed momentum loop and Berry curvature in a part of the Brillouin zone.

Figures

Equations

$${\varphi }_{\mathbf{R}m}(\mathbf{r})\equiv {\varphi }_m(\mathbf{r}-\mathbf{R})\equiv {\phi }_{i\mu }(\mathbf{r}-\mathbf{R}-{\mathbf{\tau }}_i)$$

$$H_{mn}(\mathbf{R})=⟨{\varphi }_{0m}|\hat{H}|{\varphi }_{\mathbf{R}n}⟩$$

$$H_{mn}(\mathbf{k})=\sum _{\mathbf{R}}{e}^{i\mathbf{k}\cdot \mathbf{R}}{H}_{mn}(\mathbf{R})$$

$$H_{mn}(\mathbf{k})=\sum _{\mathbf{R}}{e}^{i\mathbf{k}\cdot (\mathbf{R}+{\mathbf{\tau }}_m-{\mathbf{\tau }}_n)}{H}_{mn}(\mathbf{R})$$

$${{\mathbb{R}}^{\prime }}_1=U_{11}{\mathbb{R}}_1+U_{12}{\mathbb{R}}_2+U_{13}{\mathbb{R}}_3$$

$${{\mathbb{R}}^{\prime }}_2=U_{21}{\mathbb{R}}_1+U_{22}{\mathbb{R}}_2+U_{23}{\mathbb{R}}_3$$

$${{\mathbb{R}}^{\prime }}_3=U_{31}{\mathbb{R}}_1+U_{32}{\mathbb{R}}_2+U_{33}{\mathbb{R}}_3$$

$${{\mathbb{R}}^{\prime }}_1\cdot ({{\mathbb{R}}^{\prime }}_2\times {{\mathbb{R}}^{\prime }}_3)={\mathbb{R}}_1\cdot ({\mathbb{R}}_2\times {\mathbb{R}}_3)$$

$$H_{mn}^{slab}({\mathbb{k}}_{\parallel })=\sum _{||R||}{e}^{i{\mathbb{k}}_{\parallel }\cdot \mathbb{R}}{H}_{mn}^{slab}(\mathbb{R})$$

$$\displaystyle H^{slab}_{mn}(\mathbb{k}_{\parallel})=\left(\begin{array}[]{cccc}H_{mn,11}(\mathbb{k}_{\parallel})&H_{mn,12}(\mathbb{k}_{\parallel})&\cdots&H_{mn,1n_{s}}(\mathbb{k}_{\parallel})\\ H_{mn,21}(\mathbb{k}_{\parallel})&H_{mn,22}(\mathbb{k}_{\parallel})&\cdots&H_{mn,2n_{s}}(\mathbb{k}_{\parallel})\\ \vdots&\vdots&\ddots&\vdots\\ H_{mn,n_{s}1}(\mathbb{k}_{\parallel})&H_{mn,n_{s}2}(\mathbb{k}_{\parallel})&\cdots&H_{mn,n_{s}n_{s}}(\mathbb{k}_{\parallel})\end{array}\right)$$

$$H_{mn,ij}({\mathbb{k}}_{\parallel })=\sum _{\mathbb{R}=\left\{{{\mathbb{R}}^{\prime }}_1,{{\mathbb{R}}^{\prime }}_2,(i-j){{\mathbb{R}}^{\prime }}_3\right\}}{e}^{i{\mathbb{k}}_{\parallel }\cdot \mathbb{R}}{H}_{mn}(\mathbb{R})$$

$$|n{k}_x{l}_y⟩=\frac{1}{2\pi }{\int }_0^{2\pi }d{k}_y{e}^{-i{k}_y{l}_y}|{\psi }_{n\mathbb{k}}⟩$$

$${\bar{y}}_n(k_x)$$

$$=$$

$${\bar{y}}_n(k_x)$$

$${\bar{y}}_n(k_x)=-\frac{1}{2\pi }\text{Im}\ln {\lambda }_n$$

$${\varphi }_n={\oint }_{\mathcal{C}}{\mathbf{A}}_n\cdot d\mathbf{k}$$

$${\Omega }_{n,\alpha \beta }$$

$$=$$

$$\varphi ={\oint }_{\mathcal{C}}\mathrm{Tr}[\text{mathbfcal}A]\cdot d\mathbf{k}$$

$${\Omega }_{mn,\alpha \beta }$$

$$=$$

$${\Omega }_{\alpha \beta }=\mathrm{Tr}{\Omega }_{mn,\alpha \beta }$$

$$\varphi =-\sum _j\mathrm{Im}\ln \det M^{(j)}=-\mathrm{Im}\ln \prod _j\det M^{(j)}$$

$${\alpha }_i={\alpha }_{i-1}(\omega -{\epsilon }_{i-1}{)}^{-1}{\alpha }_{i-1}$$

$${\beta }_i={\beta }_{i-1}(\omega -{\epsilon }_{i-1}{)}^{-1}{\beta }_{i-1}$$

$${\epsilon }_i={\epsilon }_{i-1}+{\alpha }_{i-1}(\omega -{\epsilon }_{i-1}{)}^{-1}{\beta }_{i-1}+{\beta }_{i-1}(\omega -{\epsilon }_{i-1}{)}^{-1}{\alpha }_{i-1}$$

$${\epsilon }_i^s={\epsilon }_{i-1}^s+{\alpha }_{i-1}(\omega -{\epsilon }_{i-1}{)}^{-1}{\beta }_{i-1}$$

$${\tilde{\epsilon }}_i^s={\tilde{\epsilon }}_{i-1}^s+{\beta }_{i-1}(\omega -{\epsilon }_{i-1}{)}^{-1}{\alpha }_{i-1}$$

$$G_s({\mathbb{k}}_{\parallel },\omega )\simeq (\omega -{\epsilon }_n^s{)}^{-1}$$

$${\tilde{G}}_s({\mathbb{k}}_{\parallel },\omega )\simeq (\omega -{\tilde{\epsilon }}_n^s{)}^{-1}$$

$$G_b({\mathbb{k}}_{\parallel },\omega )\simeq (\omega -{\epsilon }_n{)}^{-1}$$

$$A({\mathbb{k}}_{\parallel },\omega )=-\frac{1}{\pi }\underset{\eta \to 0^{+}}{\lim }\mathrm{Im}\mathrm{Tr}{G}_s({\mathbb{k}}_{\parallel },\omega +i\eta )$$

$$\mathbb{S}({\mathbb{k}}_{\parallel },\omega )=-\frac{1}{\pi }\underset{\eta \to 0^{+}}{\lim }\mathrm{Im}\mathrm{Tr}\left[\mathbf{\sigma }{G}_s({\mathbb{k}}_{\parallel },\omega +i\eta )\right]/A({\mathbb{k}}_{\parallel },\omega )$$

Full text

PROGRAM SUMMARY

Program title: WannierTools

Catalogue identifier: TO BE DONE

Program summary URL: TO BE DONE

Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland

Licensing provisions: GNU General Public Licence 3.0

No. of lines in distributed program, including test data, etc.: 17118 lines

No. of bytes in distributed program, including test data, etc.: 7421589 bytes

Distribution format: tar.gz

Programming language: Fortran 90

Computer: any architecture with a Fortran 90 compiler

Operating system: Unix, Linux, Mac OS X

Has the code been vectorised or parallelized?: Yes

RAM: Variable, depends on the complexity of the problem

Classification: 40

External routines/libraries used:

BLAS (http://www/netlib.org/blas)

• 2.

LAPACK (http://www.netlib.org/lapack)

Nature of problem: Identifying the topological classifications of crystalline systems including insulators, semimetals, metals, and studying the electronic properties of the related slabs and ribbons systems.

Solution method: Tight-binding method is a good approximation for solid system. Based on that, Wilson loop is used for topological phase classification. The iterative Green’s function is used for obtaining the surface state spectrum.

Running time: Depends on the complexity of the problem

1 Introduction

Novel topological states have attracted much attention during the past decades. They give a lot of opportunities to explore new physics and to realize new quantum devices [2, 3]. Since quantum spin hall effect (QSHE) was predicted in 2005 [4] , and realized in HgTe/CdTe quantum wells in 2007 [5, 6], more and more topological novel states were discovered, such as 3D topological insulators [7, 8, 9], Dirac [10, 11], Weyl [12, 13, 14] semimetals, Hourglass fermions [15], Nodal line semimetals  [16, 17, 18] and nodal chain metals [19] et al… Nowadays, a lot of new topological phases are emerging, and more and more materials are identified to be topologically non-trivial.

One important feature of topological materials is topologically protected surface [20] states, which are robust against the white noise disorder and have novel transport properties [6] due to spin-momentum locked features. In experiments, the surface states and transport properties are detectable, therefore, taken as an evidence for non-trivial topology in bulk band structures  [6, 7, 21]. The logic here is that the bulk-edge correspondence principle, which tells that there are topologically protected $d-1$-dimension edge states if the topological property is non-trivial in $d$-dimension bulk. In the theoretical part, besides calculation of surface-state spectrum and transport properties, we could also calculate topological numbers such as the $\mathbb{Z}_{2}$ number [22], Chern number [23] and Wilson loops [24] (identical to Wannier charge centers [25]) to study the topology of energy band structures of a material directly.

At present, there are more and more groups joining the investigation of novel topological materials. However, there are only a few software packages that can be used so far. Z2pack [26] is a package that uses the Wannier charge centers (WCCs) to classify the topological properties of real compounds. By using Z2pack, the WCCs can be obtained either from tight-binding (TB) model or from first-principle packages, such as VASP, ABINIT and Quantum-espresso. PythTB  [27] is another software package with a Python implementation of TB models. There are many functions to build and solve TB models of the electronic structure of systems for arbitrary dimensional systems (crystals, slabs, ribbons, clusters, etc.) in PythTB, and it can compute Berry phases and related properties.

Here we introduce a TB dependent open-source software package called WannierTools, which can be used for novel topological materials investigation. Unlike Z2pack, it can be used to calculate the surface states of materials, and, being parallelized with MPI, it is faster than PythTB. It is a user friendly and efficient single program implemented in Fortran90. The only thing needed to run it is an input file, which contains some parameters describing your systems, and a TB model written in Wannier90_hr.dat format [1]. With WannierTools, topological numbers like the $\mathbb{Z}_{2}$ numbers or Wilson loops for the bulk system can be calculated in order to explore the topological properties of a material. It can also help to search for Weyl/Dirac points or nodal loop structures in the BZ of metallic systems. There are plenty of other functions, e.g. studying the electron structure properties for slab and ribbon systems, studying the Berry curvature for bulk systems, studying the Berry phase around one momentum loop in the BZ for nodal-line systems and so on.

This paper is organized as follows. In section 2, we review briefly some basic theories related to this package. In section 3, we introduce the capabilities of this package. In section 4, we introduce the installation and basic usages. In section 5, we introduce a new topological material HfPtGe in order to show you how to use WannierTools to explore a new topological phase.

2 Methods

2.1 TB method

TB method is a semi-empirical approach to study electronic structures of solid-state systems by projecting the Hamiltonian of the system onto a series of local orbitals. There are several ways to construct TB models, such as Slater-Koster method [28], maximum localized Wannier functions(MLWF) [29], and discretization of $k\cdot p$ model [30] onto a lattice. Among these methods, the MLWF method[29] is widely used by the people who are interested in real materials simulations. MLWF is implemented in Wannier90, which has many interfaces with different first-principle software packages like VASP, WIEN2k, et al.. Therefore MLWF TB models can be automatically obtained from first-principle calculations together with Wannier90 [1].

In different TB methods, the basis functions could be mutually orthogonal or non-orthogonal. However, WannierTools is only capable of dealing with the TB models with orthogonal basis functions. Fortunately, the Wannier functions (WFs) for MLWF TB method fulfill this limitation. In this section, we give some brief introductions to the general orthogonal TB methods. The details of how to construct MLWF TB models can be found in Refs. [31, 29]

Let $i$ label the atoms, $\mu$ label the orbitals, $m$ label the combination of {$i\mu$}, $\boldsymbol{R}$ label the lattice vectors in 3D crystal, and $\boldsymbol{\tau}_{i}$ label the position of atoms in a home unit cell. The local orbital for the i’th atom centered at $\boldsymbol{R}+\boldsymbol{\tau}_{i}$ can be written as

[TABLE]

The orthogonality of orbitals requires $\langle\phi_{\boldsymbol{R}m}|\phi_{\boldsymbol{R}^{\prime}n}\rangle=\delta_{{\boldsymbol{RR^{\prime}}}}\delta_{mn}$. TB parameters of the Hamiltonian that have the translational symmetry due to Bloch theorem can be calculated via

[TABLE]

Once we have the TB Hamiltonian $H_{mn}(\boldsymbol{R})$, the Hamiltonian in k space can be obtained by a Fourier transformation (FT)  [29]. There are two conventions [27] for FTs. One is

[TABLE]

The other one is

[TABLE]

It can be demonstrated that eigenvalues for these two conventions are the same, but the eigenvectors are different. The eigenvectors of the first convention Eqn.(3) are analogous to the Bloch wave functions $\psi_{n\boldsymbol{k}}(\boldsymbol{r})$. The eigenvectors of the second convention Eqn.(4) are analogous to the periodic part of the Bloch wave functions $u_{n\boldsymbol{k}}(\boldsymbol{r})=\psi_{n\boldsymbol{k}}(\boldsymbol{r})e^{-i\boldsymbol{k}\boldsymbol{r}}$, which is of great importance in Berry phase and Berry curvature or the Wannier centers calculations. Therefore, the second convention is used in WannierTools.

According to the bulk-edge correspondence, there are topologically protected surface states if the topology of bulk energy bands is non-trivial. In order to study such surface states, we have to construct a slab system which is periodic along two directions at the surface. In practice, a new unit cell is defined with lattice vectors $\mathbb{R^{\prime}}_{1,2,3}$,

[TABLE]

Where $\mathbb{R}_{1,2,3}$ are lattice vectors of the original unit cell of the bulk system, $\mathbb{R^{\prime}}_{1}$ and $\mathbb{R^{\prime}}_{2}$ are two lattice vectors in the target slab surface, $\mathbb{R^{\prime}}_{3}$ is the other lattice vector which is out of the surface and fulfills the volume fixed condition,

[TABLE]

Since the slab system is non-periodic along the $\mathbb{R^{\prime}}_{3}$ direction, the Hamiltonian of a slab system with a 2D momentum $k_{\parallel}$ can be calculated by the following FT

[TABLE]

Where $\mathbb{R}=a^{\prime}\mathbb{R^{\prime}}_{3}+b^{\prime}\mathbb{R^{\prime}}_{2}+c^{\prime}\mathbb{R^{\prime}}_{3}$ and $||\mathbb{R}||$ is a restriction that the summation is only carried on $a^{\prime}$ and $b^{\prime}$ with different $c^{\prime}$. We label the layer index along $\mathbb{R^{\prime}}_{3}$ as $i,j$. As a consequence, the Hamiltonian of a slab system with $n_{s}$ layers can be written in the layer index matrix form,

[TABLE]

where the diagonal elements of the Hamiltonian are the intra-plane ones, and the off diagonal elements are the inter-plane ones. The element in Eq. (14) can be read explicitly as,

[TABLE]

Finally the energy band of the slab system can be obtained straightforwardly to diagonalize the Eq. (14).

By the way, the Hamiltonian for the ribbon system can be obtained in the same way as the slab system does. The difference is that there are two confined directions $\mathbb{R^{\prime}}_{1},\mathbb{R^{\prime}}_{2}$ in ribbon systems, which enlarges the size of Hamiltonian.

2.2 Wannier charge center calculation

$\mathbb{Z}_{2}$ topological number [22] and Chern number [23] are applied to classify topological properties for time-reversal invariant and time-reversal symmetry breaking systems respectively. In inversion symmetric invariant system, the $\mathbb{Z}_{2}$ topological number can be calculated by multiplying the parities of the occupied bands at time reversal invariant momenta (TRIMs) in the Brillouin zone [32]. There are several methods [33, 24, 25] to calculate the $\mathbb{Z}_{2}$ number in inversion symmetry breaking systems. Among them it was demonstrated that the Wilson loop [24] and Wannier charge centers (WCCs) [25] method are equivalent to each other, and are also valid for time-reversal symmetry breaking systems. In WannierTools, we take the algorithm presented in Refs. [25, 26]. The hybrid WFs [34] are defined as

[TABLE]

where $|\psi_{n{\mathbb{k}}}\rangle$ is the Bloch wave function. The hybrid Wannier centers are defined as

[TABLE]

where $|u_{n,k_{x},k_{y}}\rangle$ is the periodic part of Bloch function $|\psi_{n{\mathbb{k}}}\rangle$. In practice , the integration over $k_{y}$ is transformed by a summation over the discretized $k_{y}$. Eqn.(18) can be reformulated using the discretized Berry phase formula [35],

[TABLE]

where the gauge-dependent overlap matrix $M_{mn}^{(j)}=\langle u_{m,k_{x},k_{y_{j}}}|u_{n,k_{x},k_{y_{j+1}}}\rangle$ is introduced. However, the summation of the hybrid Wannier centers $\bar{y}_{n}(k_{x})$ over $k_{x}$ is gauge invariant [36], which is the total electronic polarization. As shown in Ref [35, 37], there is another way to obtain $\bar{y}_{n}(k_{x})$. Firstly, we get the ”unitary part” $\tilde{M}^{(j)}$ of each overlap matrix $M_{mn}^{(j)}$ by carrying out the single value decomposition $M=V\Sigma W^{{\dagger}}$, where V and W are unitary and $\Sigma$ is real-positive and diagonal. Then we set $\tilde{M}^{(j)}=VW^{{\dagger}}$. The eigenvalues $\lambda_{n}$ of matrix $\Lambda=\prod_{j}\tilde{M}^{(j)}$ are all of unit modulus. The hybrid Wannier centers are defined with the phases of $\lambda_{n}$ eventually,

[TABLE]

We can get the topological properties of $k_{x}-k_{y}$ plane from the evolution of $\bar{y}_{n}(k_{x})$ along a $k_{x}$ string. The detail of such classification of WCCs or Wilson loop are discussed in Refs [24, 37, 25]. More information can be found in Ref  [26].

2.3 Berry phase and Berry curvature

In this section, we give the basic formalism for computing Berry phase [38, 39] and Berry curvature [40, 41] of Bloch states. Firstly, we introduce the single band case, where the energy bands are isolated to each other. A Berry phase $\phi_{n}$ is a geometric phase associated with the phase evolution of the n’th state over a closed curve $\mathcal{C}$ in external parameter space k, defined as

[TABLE]

Where the Berry connection is $A_{n,\alpha}=i\langle u_{n,{\boldsymbol{k}}}|\partial_{\alpha}u_{n,{\boldsymbol{k}}}\rangle,\alpha=k_{x},k_{y},k_{z}$, and Berry curvature are introduced

[TABLE]

Secondly, for multi-band case, it is often to treat the occupied $N_{occ}$ bands as a joint band manifold, which is referred to as the ”non-Abelian” case. Generalizations for the formalism of Berry phase and Berry curvature from single band to multi-band case are as follows,

[TABLE]

Where the Berry connection  [39] is $\mathcal{A}_{mn,\alpha}=i\langle u_{m,{\boldsymbol{k}}}|\partial_{\alpha}u_{n,{\boldsymbol{k}}}\rangle,\alpha=k_{x},k_{y},k_{z}$, and Berry curvature is

[TABLE]

and define

[TABLE]

where $\operatorname{Tr}$ denotes a trace over the occupied bands.

In practice, the integration of Eqn.(25) is implemented on a discrete k mesh. The loop $\mathcal{C}$ is discretized into a series of closely space points $k_{j}$. Accordingly, the Berry phase becomes

[TABLE]

where the overlap matrix $M_{mn}^{(j)}$ is the same as in Eqn.(19), i.e., $M_{mn}^{(j)}=\langle u_{m,{\boldsymbol{k}}_{j}}|u_{n,{\boldsymbol{k}}_{j+1}}\rangle$.

2.4 Calculation of surface states

Theoretically, we have two methods to get surface spectrum corresponding to the bulk topology. One is that we calculate the band structure of a slab system, which was introduced in section 2.1. The other one is to calculate the surface Green’s function (SGF) for a semi-infinite system that will be introduced in this section. In the 1970s, one of the most popular GF approachs was based on the ”effective field” and transfer matrix [42, 43, 44], which are relatively of slow convergence especially near singularities. Now, the extensively used schemes to obtain the SGFs is the iterative Green’s function method developed in the 1980s [45, 46]. With an effective concept of principle layers (The layer that is large enough to guarantee that hoppings between the next nearest layers are negligible.), the iterative procedure can save quite an amount of computational time. The method [45, 46] involves replacing the principle layer by an effective two principle layers, and these effective layers interact through energy-dependent residual interactions which are weaker than those of the original ones. This replacement can be repeated iteratively until the residual interactions between the effective layers become as small as desired. Each new iteration doubles the number of the original layers included in the new effective layer. That is, after $n$ iterations, one has a chain of lattice constant $2^{n}$ times the original one, and each effective layer replacing $2^{n}$ original layers. The details of the algorithm are presented in Ref [47]. For the integrity of the paper, we list the main iterations that we reused in WannierTools. The most important parameters for iteration $i$ are the following

[TABLE]

with the initialization $\varepsilon_{0}=\tilde{\varepsilon}_{0}=\tilde{\varepsilon}_{0}^{s}=H_{00}({\mathbb{k}}_{\parallel})$, $\alpha_{0}=H_{01}({\mathbb{k}}_{\parallel})$, $\beta_{0}=H_{01}^{{\dagger}}({\mathbb{k}}_{\parallel})$, where $H_{00}({\mathbb{k}}_{\parallel})$ is the intra-hopping parameters in the principle layers, $H_{01}({\mathbb{k}}_{\parallel})$ is the inter-hopping parameters between the nearest neighbor of principle layers.

Iteration of Eqn.(31) should be converged until ${\varepsilon}_{n}\simeq{\varepsilon}_{n-1}$ and $\tilde{\varepsilon}_{n}^{s}\simeq\tilde{\varepsilon}_{n-1}^{s}$. The SGFs $G_{s}({\mathbb{k}}_{\parallel},\omega)$ and the bulk GF $G_{b}({\mathbb{k}}_{\parallel},\omega)$ can be obtained as

[TABLE]

Where $\tilde{G}_{s}$ is the SGF of the dual surface. The surface spectrum function $A({\mathbb{k}}_{\parallel},\omega)$ can be obtained from the imaginary part of SGF

[TABLE]

The spin texture of surface states can also be obtained with [48]

[TABLE]

where ${\boldsymbol{\sigma}}$ are the Pauli matrixes.

2.5 Algorithm for searching nodal points/lines

Nodal point is a gapless point between the highest valence band and the lowest conduction band. In terms of degeneracy character, nodal point could be classified into Weyl, Triple, Dirac, hyper-Dirac point with 2-fold degeneracy, 3-fold degeneracy, 4-fold degeneracy or higher degeneracy respectively. They also can be sorted into nodal point and nodal line by the connectivity between them. Searching the Weyl/Dirac points and the nodal-line structures is very important for such nodal systems. Some symmetry protected nodal points or nodal lines which are located in high symmetry lines or mirror planes are easy to find. While the other nodes which are located anywhere in the BZ need more efforts to be found. Here we introduce an algorithm trying to find all the nodal points.

Basically, node points are local minima of the energy gap function in 3D BZ. Local minima can be obtained by using some well known multidimensional minimization methods, e.g., Nelder and Mead’s Downhill Simplex Method [49], Conjugate Gradient Methods [50], Quasi-Newton Methods [51] et al.. However, the local minimum obtained from those methods depends on a initial point. One initial point gives only one local minimum. So, in order to find all the nodes, we have to choose different initial points in the whole 3D BZ. WannierTools takes a uniform mesh of the 1st BZ as a set of initial points for the Nelder and Mead’s Downhill Simplex Method [49]. Eventually, the nodes will be selected out from a set of local minima. It is easy to check the convergence of the number of nodes by increasing the initial point mesh. This algorithm is very suitable for high throughput search of new Weyl, Dirac semimetals and nodal-line metals. It has been checked to be very efficient to find Weyl points in WTe2 [14, 52], MoTe2 [53].

3 Capabilities of WannierTools

There are two kinds of tasks that WannierTools can do to study novel topological materials. a. One is to get the topology of materials’ band structure. b. The other one is to explore the properties of surface states corresponding to the bulk topology. For part a, we need to study the bulk band structure, 3D Fermi surface, density of state (DOS) to check whether the bulk material is a band insulator or a metal. Further, WCCs calculations are applied to get the $\mathbb{Z}_{2}$ topological index or Chern number for band insulators, the nodes searching algorithm and the energy gap function calculation are applied to search for Weyl/Dirac point positions or nodal-line structures; The Berry phase and Berry curvature calculations are also aided to the classification of topology. After the topological classification is done, one can turn to part b, which means to study the bulk topology related properties, such as joint density of state (JDOS), which is related to the optical conductivity [54], electronic structure of the slab and wire systems, spin-texture of the surface states, Quasi-particle interference (QPI) pattern of surface states et al.. These two main capabilities of WannierTools are listed in Table 3 and Table 3. The meaning of control flags in Table 3 and Table 3 is illustrated in the documentation, which is distributed with WannierTools.

[htp] Main capabilities of WannierTools: Bulk topology studies

Control flag in wt.in Description

BulkBand_calc Energy bands for a 3D bulk system.

BulkFS_calc 3D Fermi surface in 1st BZ.

FindNodes_calc Locate Weyl, Dirac point positions and nodal line structures in 1st BZ.

BulkGap_plane_calc Gap function in a 3D k plane.

Wanniercenter_calc WCCs [25] for a 3D k plane.

BerryPhase_calc Berry phase for a closed path in 3D k space.

BerryCurvature_calc Berry curvature in a 3D k plane.

[htp] Main capabilities of WannierTools: related responses from the bulk topology

Control flag in wt.in Description

Dos_calc Density of state of a 3D bulk system

JDos_calc Joint density of state [54] of a 3D bulk system

SlabBand_calc Energy bands of a 2D slab system.

WireBand_calc Energy bands of a 1D ribbon system.

SlabSS_calc Surface state spectrum along some k lines at different energies.

SlabArc_calc Surface state spectrum in the 2D BZ at a fixed energy.

SlabQPI_calc Quasi-particle interference (QPI) [55] pattern of surface state.

SlabSpintexture_calc Spin texture [56] of surface state.

4 Installation and usage

In this section, we will show how to install and use the WannierTools software package.

4.1 Get WannierTools

WannierTools is an open source free software package. It is released on Github under the GNU General Public Licence 3.0 (GPL), and it can be downloaded directly from the public code repository:

https://github.com/quanshengwu/wannier_tools.

4.2 Build WannierTools

To build and install WannierTools, a Fortran 90 compiler, BLAS, and LAPACK linear algebra libraries are needed. An MPI-enabled Fortran 90 compiler is needed if you want to compile a parallel version. WannierTools can be successfully compiled using the state-of-art Intel Fortran compiler. Most of the MPI implementations, such as MPICH, MVAPICH and Intel MPI are compatible with WannierTools. The downloaded WannierTools software package is likely a compressed file with a zip or tar.gz suffix. One should uncompress it firstly, then move into the wannier_tools/soc folder and edit the Makefile file to configure the compiling environment. It is noteworthy that one should set up the Fortran compiler, BLAS and LAPACK libraries manually by modifying the following lines in Makefile file according to the user’s particular system.

f90 = libs =

Once the compiling environment is configured, the executable binary wt.x will be compiled by typing the following command in the current directory (wannier_tools/soc)

$ make

4.3 Running WannierTools

Before running WannierTools, the user must provide two files wannier90_hr.dat111The name of this TB file will be shown in the wt.in file. So it could be any name you like. and wt.in. The file called wannier90_hr.dat containing the TB parameters has fixed format which is defined in software Wannier90 [1]. It can be generated by the software Wannier90 [1], or generated by users with a toy TB model, or generated from a discretization of $k\cdot p$ model onto a cubic lattice. The other file wt.in is the master input file for WannierTools. It is designed to be simple and user friendly. The details of wt.in are described comprehensively in the documentation that contained within the WannierTools distribution. An example file is provided in A.

After putting wt.in and wannier90_hr.dat in the same folder, one can run it in single processor in the same folder like this

$ wt.x &

or in multiprocessor

$ mpirun -np 4 wt.x &

Some important information during the running process are written in WT.out, from which, you can check the running status. After the whole program is done, you would obtain two kinds of files other than WT.out. One is the data file suffixed with dat. and the other one is a plotting script for software gnuplot suffixed with gnu. You can get nice plots with gnuplot [57]. Taking a bulk band structure calculation from the examples, two files bulkek.dat and bulkek.gnu are accomplished after a successful running of WannierTools. A band structure plot bulkek.png will be generated with the following command

gnuplot bulkek.gnu

5 Examples

In the past few years, WannierTools has been successfully applied in many projects, such as finding type-II Weyl semimetals WTe2 [14, 52], MoTe2 [53], triple point metals [58] ZrTe, TaN, nodal chain metals IrF4 [19], topological phase in InAs1-xSbx  [59] et al.. Besides, more and more groups notice this package, and becoming users. There are several examples in the wannier_tools/examples directory like Bi2Se3, WTe2, IrF4. MLWF TB Hamiltonians for those materials and the necessary input files for generating those Hamiltonians can be downloaded from the Github repository [60]. The detailed hands-on tutorials for those examples are listed in the wiki of Github  [61]. In this paper, a new series of topological materials called ternary silicides and ternary germanides TiPtSi, ZrPtSi, ZrPtGe, HfPtSi and HfPtGe  [62, 63] are exhibited as an example to show how to study topological properties of new materials with WannierTools.

5.1 Crystal structure and Band structure

Ternary silicides and ternary germanides are crystalized with the orthorhombic TiNiSi type structure in a nonsymmorphic orthorhombic space group No.62 (Pnma) [64, 62], containing three glide reflections $G_{x}=\{m_{x}|\frac{1}{2},\frac{1}{2},\frac{1}{2}\}$, $G_{z}=\{m_{z}|\frac{1}{2},0,\frac{1}{2}\}$, $\tilde{G}_{y}=\{m_{y}|0,\frac{1}{2},0\}$  [65], three screw rotations $S_{x}=\{R_{2x}|\frac{1}{2},\frac{1}{2},\frac{1}{2}\}$, $S_{z}=\{R_{2z}|\frac{1}{2},0,\frac{1}{2}\}$, $S_{y}=\{R_{2y}|0,\frac{1}{2},0\}$ and an inversion symmetry $I$. These materials seem to be interesting systems in the search for new superconducting intermetallic compounds [63]. Due to the same crystal structure and similar chemical properties, these materials show very similar band structures. In this paper, we are only focusing on HfPtGe compound.

As mentioned in section 4.3, we need a TB model of HfPtGe for WannierTools. Firstly, a first-principle calculation within Vienna Ab initio Simulation Package (VASP)  [66, 67] using Gamma centered K-points mesh $8\times 13\times 7$ and energy cut 360eV for plane wave expansions was performed. Then the band structure and the partial density of states (PDOS) shown in Fig 1c and Fig 1d were analyzed, where the PDOS indicate that the relevant orbitals close to the Fermi level are dominated by Hf 5d orbitals and Pt 5d orbitals, besides, they are also hybridized with Ge 4p and Pt s orbitals. In the end, a 112-band MLWF TB Hamiltonian with Hf 5d, Pt 6s 5d and Ge 4p as projectors are constructed with Wannier90. Fig.1c shows that the band structures calculated from the MLWF TB model are quite fitted to the first-principle calculated band structures. After the successful construction of MLWF TB model, the master input file wt.in is needed, which is attached in the A.

5.2 Energy gap shape

There are two ”Dirac like” cones along Y-$\Gamma$-Z direction which are shown in the band structure plot Fig. 1 (c). From the literature, we know that such cones could be originated from a nodal line structure without SOC [17]. Calculations of the energy gap at $\text{k}_{\text{y}}=0$ plane both for without SOC and with SOC cases were performed to study the details of positions of these ”Dirac like” cones. The results are shown in Fig.2a and Fig.2b, where there is a gapless nodal line protected by the $\tilde{G_{y}}$ mirror symmetry [17] in the $\text{k}_{\text{y}}=0$ plane without SOC , while, the gapless line will be gapped when SOC turns on. One thing should mentioned is that the SOC strength of this material is very weak, and the opened gap is not very big. The smallest gap is about 0.1meV where the k points are located at ($\pm$0.4, 0.0, $\pm$0.229).

We can further study the nodal line distribution in momentum and energy space, which shown in Fig.2c. It is clearly shown that the nodal line is not in the same energy plane. There are six nodes crossing the Fermi level. In such case, there are electron pockets and hole pockets which link together on the Fermi surface at the same time, which are shown in Fig 2 (d). Those compensated hole and electron pockets will cause extremely large positive magnetoresistance [68].

The related settings in the master input file wt.in for this section are as follows

&CONTROL BulkGap_plane_calc = T BulkFS_calc = T /

&PARAMETERS
NK1= 101
NK2= 101
NK3= 101
/

KPLANE_BULK
-0.50  0.00 -0.50   ! Original point for 3D k plane
 1.00  0.00  0.00   ! The first vector to define 3d k space plane
 0.00  0.00  1.00   ! The second vector to define 3d k space plane

5.3 Wannier charge center

From the gap shape calculation for the whole BZ in the SOC case, we can conclude that HfPtGe is a semimetal with a continuous finite energy gap between electron-like and hole-like bands. Similar to classification of band insulators [69], $\mathbb{Z}_{2}$ topological indices $(\nu_{0},\nu_{1}\nu_{2}\nu_{3})$ are still appropriate for such a semimetal. $\mathbb{Z}_{2}$ number of a bulk material can be obtained through calculations of WCCs in six time reversal invariant planes $\text{k}_{\text{x}}=0,\pi$, $\text{k}_{\text{y}}=0,\pi$ and $\text{k}_{\text{z}}=0,\pi$ plane. The results calculated by WannierTools are shown in Fig. 3. It shows that the $Z_{2}$ invariant numbers are 1 for $\text{k}_{\text{x}}=0,\text{k}_{\text{y}}=0,\text{k}_{\text{z}}=0$ plane, while zeros for other planes. Eventually, the topological index is $(1,000)$, which indicates that HfPtGe is a ”strong” topological material in all three reciprocal lattice directions. The related settings for this section in wt.in are as follows.

&CONTROL
Z2_3D_calc = T
/

5.4 Surface state spectrums

Due to the bulk-edge correspondence, there should be topologically protected surface states of any cuts of surface for a strong topological material [69]. Here we study the (010) surface which is shown as a grey plane in Fig.1a. The surface state spectrums calculated by WannierTools are shown in Fig. 4a,4b. For a 3D strong topological insulator, there is surface Dirac cone at the $\Gamma$ point [8]. Indeed, there is a Dirac like cone of HfPtGe at the $\Gamma$ point. However, the dispersion of this cone is highly anisotropic and even tiled in the momentum space. Such tiled cone is result from that the nodal line is located at different energies. Fig. 4b shows a $\text{E}=\text{E}_{\text{F}}$ iso-energy plot of the surface state spectrum. From this plot, we can learn that the surface states originate from the k points that have the smallest gap. The related settings for this section in wt.in are as follows.

&CONTROL
SlabSS_calc  = T
SlabArc_calc = T
/

&PARAMETERS
E_arc = 0.0    ! Fixed energy Fermi arc calculation
/

KPATH_SLAB
4        ! numker of k lines for a slab system
X -0.50  0.00  G  0.00 0.00  ! k path for a slab system
G  0.00  0.00  Z  0.00 0.50
Z  0.00  0.50  R -0.50 0.50
R -0.50  0.50  X -0.50 0.00

KPLANE_SLAB
-0.5 -0.5
 1.0  0.0
 0.0  1.0

5.5 Surface state spin texture

The key ingredient to generate topological non-trivial properties is the SOC interaction, which have a similar behavior as a Lorentz force in quantum Hall states. The SOC interaction will make a spin and a momentum locked to each other, forming a spin texture in momentum space. For different topological phase, the spin texture will be different [70]. In topological insulator Bi2Se3, the spin texture is Dirac type. The spin texture calculated by WannierTools for HfPtGe is shown in Fig.4c. The related settings for spin texture calculations in wt.in are as follows.

&CONTROL
SlabSpintexture_calc = T
/

KPLANE_SLAB
-0.50 -0.50
 1.00  0.00
 0.00  1.00

6 Conclusions

In conclusion, we presented an open-source software package called WannierTools. It is very user-friendly and is written in Fortran90, using MPI techniques to get excellent performance in computer cluster. We showed how to use this software package to identify the topological properties for a new material and to get the surface state spectrum which can be compared with experimental data. As an example, we explored a new topological material HfPtGe, which was identified as a Dirac nodal line semimetal.

Acknowledgments

QSW would like to thank Xi Dai, Zhong Fang, Lei Wang, Li Huang, Dominik Gresch, Zhida Song for helpful discussions. Especially, QSW appreciates Rui Yu and Haijun Zhang for their kindly helps at the beginning of this project. QSW, AAS, MT were supported by Microsoft Research, and the Swiss National Science Foundation through the National Competence Centers in Research MARVEL and QSIT, ERC Advanced Grant SIMCOFE. QSW was also supported by the National Natural Science Foundation of China (11404024). SNZ was supported by NSF-China under Grants No. 11074174. HFS was supported by National High Technology Research and Development Program of China under Grant 2015AA01A304, and Science Challenge Project No. JCKY2016212A502. This job was started in IOP CAS, finished in ETH Zurich.

Note added.—Recently, Ref.  [71] appeared, discussing topological properties of the same class of HfGePt.

Appendix A wt.in for HfPtGe

&TB_FILE Hrfile = "wannier90_hr.dat" /

&CONTROL BulkBand_calc = T BulkFS_calc = T BulkGap_plane_calc = T Z2_3D_calc = T SlabSS_calc = T SlabArc_calc = T SlabSpintexture_calc = T /

&SYSTEM NumOccupied = 64 SOC = 1 E_FERMI = 8.4551 /

&PARAMETERS E_arc = 0.0 OmegaNum = 100 OmegaMin = -0.6 OmegaMax = 0.6 Nk1 = 101 Nk2 = 201 Nk3 = 101 NP = 2 Gap_threshold = 0.05 /

LATTICE Angstrom 6.6030000 0.0000000 0.0000000 0.0000000 3.9500000 0.0000000 0.0000000 0.0000000 7.6170000

ATOM_POSITIONS 12 Direct Hf 0.029900 0.250000 0.186300 Hf 0.470100 -0.250000 0.686300 Hf -0.029900 0.750000 -0.186300 Hf 0.529900 0.250000 0.313700 Ge 0.750600 0.250000 0.621500 Ge -0.250600 -0.250000 1.121500 Ge -0.750600 0.750000 -0.621500 Ge 1.250600 0.250000 -0.121500 Pt 0.142500 0.250000 0.561700 Pt 0.357500 -0.250000 1.061700 Pt -0.142500 0.750000 -0.561700 Pt 0.642500 0.250000 -0.061700

PROJECTORS 5 5 5 5 3 3 3 3 6 6 6 6 Hf dxy dyz dxz dx2-y2 dz2 Hf dxy dyz dxz dx2-y2 dz2 Hf dxy dyz dxz dx2-y2 dz2 Hf dxy dyz dxz dx2-y2 dz2 Ge px py pz Ge px py pz Ge px py pz Ge px py pz Pt s dxy dyz dxz dx2-y2 dz2 Pt s dxy dyz dxz dx2-y2 dz2 Pt s dxy dyz dxz dx2-y2 dz2 Pt s dxy dyz dxz dx2-y2 dz2

SURFACE 1 0 0 0 0 1 0 -1 0

KPATH_BULK 9 Y -0.50000 0.00000 0.00000 G 0.00000 0.00000 0.00000 G 0.00000 0.00000 0.00000 Z 0.00000 0.00000 0.50000 Z 0.00000 0.00000 0.50000 T -0.50000 0.00000 0.50000 T -0.50000 0.00000 0.50000 Y -0.50000 0.00000 0.00000 Y -0.50000 0.00000 0.00000 S -0.50000 0.50000 0.00000 S -0.50000 0.50000 0.00000 X 0.00000 0.50000 0.00000 X 0.00000 0.50000 0.00000 U 0.00000 0.50000 0.50000 U 0.00000 0.50000 0.50000 R -0.50000 0.50000 0.50000 R -0.50000 0.50000 0.50000 S -0.50000 0.50000 0.00000

KPATH_SLAB 4 X -0.50 0.00 G 0.00 0.00 G 0.00 0.00 Z 0.00 0.50 Z 0.00 0.50 R -0.50 0.50 R -0.50 0.50 X -0.50 0.00

KPLANE_SLAB -0.5 -0.5 1.0 0.0 0.0 1.0

KPLANE_BULK -0.50 0.00 -0.50 1.00 0.00 0.00 0.00 0.00 1.00

KCUBE_BULK -0.50 -0.50 -0.50 1.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 1.00