Z2Pack postprocessor¶
The Z2Pack interface tutorial¶
This tutorial aims at showing how to use the Z2Pack application with an interface to ABINIT, to compute the Z2 topological invariant and find topologically nontrivial insulators.
You will learn how to use Wannier90 and Z2Pack to follow Wannier charge centers which allow to compute the Z2 topological invariant. You will be reproducing some of the results obtained in [BrousseauCouture2020], where it it shown that BiTeI is a trivial insulator at 0GPa but undergoes a transition to a nontrivial insulating phase under isotropic pressure. At 5GPa, the material is a nontrivial topological insulator.
You already should know about the interface between ABINIT and Wannier90 (see the tutorial Wannier90). This tutorial should take about 1 hour. It is important to note that the examples in this tutorial are not converged, they are just examples to show how to use the code.
Note
Supposing you made your own installation of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit toplevel directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.
In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:
export ABI_HOME=Replace_with_absolute_path_to_abinit_top_level_dir # Change this line
export PATH=$ABI_HOME/src/98_main/:$PATH # Do not change this line: path to executable
export ABI_TESTS=$ABI_HOME/tests/ # Do not change this line: path to tests dir
export ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/ # Do not change this line: path to pseudos dir
Examples in this tutorial use these shell variables: copy and paste
the code snippets into the terminal (remember to set ABI_HOME first!) or, alternatively,
source the set_abienv.sh
script located in the ~abinit directory:
source ~abinit/set_abienv.sh
The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.
To execute the tutorials, create a working directory (Work*
) and
copy there the input files of the lesson.
Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel with MPI, see the topic on parallelism.
1. Summary of Z2Pack in ABINIT¶
Z2Pack is a code that enables the calculation of topological invariants from model Hamiltonians, as well as from firstprinciples calculations. It is highly recommended to read the following papers to understand its basics: [Gresch2017] and [Gresch2018].
Z2Pack starts from a oneparticle Hamiltonian H(\textbf{k}) which describes the electronic states in a crystalline material. Its eigenvectors are socalled Bloch states  \psi_{n,\textbf{k}} \rangle. They are superpositions of plane waves with wave vector \textbf{k}. Bloch states can be written as  \psi_{n, \textbf{k}} \rangle = e^{i\textbf{k} \cdot \textbf{r}}  u_{n, \textbf{k}} \rangle where u_{n, \textbf{k}} \rangle is cellperiodic. Taking into account the shape of the Bloch states leads to a topological classification of materials.
The topological classification of materials is based on a topological invariant called the Chern number. For a closed, orientable twodimensional surface S in reciprocal space, the Chern number invariant C can be defined in terms of the cellperiodic states in a set B of bands as C = \frac{i}{2\pi} \sum_{n\in B} \int_S \nabla_{\textbf{k}} \wedge \langle u_{n, \textbf{k}}  \nabla_{\textbf{k}}  u_{n, \textbf{k}} \rangle \cdot d\textbf{S}.
There are various ways to calculate the Chern number that are easier numerically. One involves dividing the surface S into small segments, which correspond to calculating the Berry phase for closed loops in the Brillouin zone. These can then be expressed as Wilson loops, for which the eigenvalues are connected to the total Berry phase.
Another approach used by Z2Pack is to compute the socalled Wannier charge centers, which is based on the notion of Wannier orbitals. These are given by Fourier transforming the Bloch states \textbf{R}n\rangle = \frac{V}{(2\pi)^d} \int_{BZ} e^{i\textbf{k} \cdot \textbf{R}}  \psi_{n, \textbf{k}} \rangle d\textbf{k} where d is the dimensionality of the system and V the unit cell volume. These orbitals can be changed by a Gauge transformation which affects their localization and position in real space. To compute topological invariants, hybrid Wannier orbitals are introduced: they are Fourier transforms performed only in one spatial direction, for example \textbf{R}_x, k_y, k_z; n \rangle = \frac{a_x}{2\pi} \int^{\pi/a_x}_{\pi/ax} e^{ik_xR_x}  \psi_{n \textbf{k}} \rangle.
The average position of such orbital depend on the remaining reciprocal space variables: \bar{x}_n(k_y, k_z) = \langle 0, k_y, k_z; n  \hat{x}  0, k_y, k_z; n\rangle. This quantity is called the hybrid Wannier charge center and it is directly related to the Berry phase as C = \frac{1}{a} \sum_n \bar{x}_n. This procedure of segmenting the surface is illustrated in [Gresch2018] Figure 6.
The equivalence between hybrid Wannier charge centers and the Berry phase gives rise to a physical interpretation of the Chern number as a charge pumping process, where each cycle of k_x moves the charge by C unit cells.
For timereversal invariant systems, the total Chern number is zero as the Chern numbers of Kramers partners cancel each other. We will rather compute the Z2 invariant, which can only take values of 0 or 1. Z2Pack can handle this calculation in the following way. By fixing one reciprocal space coordinate at a timereversal invariant value (i.e. k_i=0 or k_i=0.5 in reduced coordinates), we can define a surface in reciprocal space which contains its own timereversal image. Using the hybrid Wannier charge centers computed on this surface, we can then define a surface invariant \Delta\in\{0,1\}, which can be intuitively understood as the Chern number mod 2 of one subset of Kramers partners. More details about how Z2Pack handles this calculation will be provided later.
To compute the surface invariant \Delta, Z2pack needs 2 inputs: a description of the system and the surface on which the invariant should be calculated. In the interface with ABINIT, Z2Pack will start from a computed ground state and perform the Wannierization on lines in the Brillouin zone. At the end, it can sum mod 2 the surface invariants obtained on the two timereversal invariant surfaces to evaluate the Z2 invariant.
Before we actually compute the Z2 invariant, we will take a look at another indication of nontrivial topology: band inversion.
2. Pressure dependent band inversion in BiTeI¶
Before beginning, you might consider working in a different subdirectory as for the other tutorials. Follow these instructions which take you to the right folder, create a subfolder called Work_z2pack and copy the first calculation:
cd $ABI_TESTS/tutoplugs/Input
mkdir Work_z2pack && cd Work_z2pack
cp ../tz2_2.abi .
You are now ready to run abinit by using (here for four cores  other number of cores or sequential runs are possible as well) :
mpirun n 4 abinit tz2_2.abi > log 2> err &
Let us examine the input file tz2_2.abi, while the calculation is running:.
# Ground state calculation for BiTeI at 0 GPa autoparal 1 pp_dirpath "$ABI_PSPDIR/Pseudodojo_nc_fr_04_pbe_standard_psp8/" pseudos "Bi.psp8,Te.psp8,I.psp8" ndtset 6 udtset 2 3 # Unit cell ntypat 3 znucl 83 52 53 natom 3 typat 1 2 3 acell1? 8.3956 8.3956 13.7627 rprim1? 0.86602540378 0.5 0.0 0.86602540378 0.5 0.0 0.0 0.0 1.0 xred1? 0.0000000000E+00 0.0000000000E+00 1.953E03 # Bi 6.6666666667E01 3.3333333333E01 2.36867E01 # Te 3.3333333333E01 6.6666666667E01 7.11179E01 # I acell2? 8.0007E+00 8.0007E+00 1.23580E+01 rprim2? 0.86602540378 0.5 0.0 0.86602540378 0.5 0.0 0.0 0.0 1.0 xred2? 0.0000000000E+00 0.0000000000E+00 4.186E03 # Bi 6.6666666667E01 3.3333333333E01 2.72237E01 # Te 3.3333333333E01 6.6666666667E01 6.73578E01 # I # Ground state calculation ecut 15 # unconverged nspinor 2 # WFK are spinors kptopt?1 1 ngkpt?1 6 6 6 nshiftk?1 1 shiftk?1 0.0 0.0 0.5 prtden?1 1 tolvrs?1 1.0e10 # unconverged nstep?1 30 nband?1 40 nbdbuf 2 # Buffer to ease convergence prtwf 0 diemac 12 enunit 1 # Full bandstructure iscf?2 2 getden?2 1 kptopt?2 7 ndivsm?2 5 kptbounds?2 0.5 0.0 0.0 #M 0.0 0.0 0.0 #Gamma 1/3 1/3 0.0 #K 1/3 1/3 0.5 #H 0.0 0.0 0.5 #A 0.5 0.0 0.5 #L 1.0 0.0 0.0 #Gamma 1.0 0.0 0.5 #A tolwfr?2 1.0e6 # unconverged nband?2 50 # Bandstructure and fatbands on HAL path iscf?3 2 getden?3 2 kptopt?3 2 prtdos?3 3 ndivsm?3 15 kptbounds?3 1/6 1/6 0.5 #0.5H 0.0 0.0 0.5 #A 0.25 0.0 0.5 #0.5L tolwfr?3 1.0e6 # unconverged nband?3 50 ratsph?3 2.5 1.9 2.15 ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [paral_info] #%% max_nprocs = 10 #%% nprocs_to_test = 10 #%% [NCPU_10] #%% files_to_test = #%% tz2_2_MPI10.abo, tolnlines=5, tolabs=7.00e05, tolrel=1.10e+00, fld_options=medium #%% [extra_info] #%% authors = O. Gingras, V. BrousseauCouture #%% keywords = wannier90, z2pack, fatbands #%% description = Test interface with Z2Pack (NC pseudopotentials) #%%<END TEST_INFO> #
This input file should look familiar, as it involves features of ABINIT that you probably already used. It is constructed into two loops of three datasets. The first big loop uses the system at 0 GPa while the second is performed at 5 GPa. In each of these big loops, the first computation calculates the ground state on a coarse grid, the second one computes the band structure along highsymmetry paths in the entire Brillouin zone and the third one is the calculation of fatbands which outputs the orbital character of the different bands, only along a specific path.
Let us compare the band structure. We use a small python script based on abipy:
#!/usr/bin/env python from abipy.abilab import abiopen import abipy.data as abidata import matplotlib.pyplot as plt import numpy as np flist = ["tz2_2o_DS12_GSR.nc", "tz2_2o_DS22_GSR.nc"] flist_gap = ["tz2_2o_DS13_FATBANDS.nc", "tz2_2o_DS23_FATBANDS.nc"] gspec = dict(width_ratios=[2, 1], height_ratios=[1, 1]) fig, ax = plt.subplot_mosaic([['full0', 'gap0'], ['full5', 'gap5']], gridspec_kw=gspec, figsize=(12, 8), constrained_layout=True) ax_list = ['full0', 'full5'] ax_gap_list = ['gap0', 'gap5'] for j in range(2): full = abiopen(flist[j]).ebands full.set_fermie_to_vbm() full.plot(with_gaps=False, ylims=(3.0, 3.0), ax=ax[ax_list[j]], show=False) gap = abiopen(flist_gap[j]).ebands gap.set_fermie_to_vbm() gap.plot(with_gaps=True, ylims=(1.0, 1.0), ax=ax[ax_gap_list[j]], show=False) ax['full5'].set_title('5 GPa') ax['full0'].set_title('0 GPa') ax['gap5'].set_xticks([0, 17, len(gap.kpoints)1]) ax['gap5'].set_xticklabels(['0.5H', 'A', '0.5L']) ax['gap0'].set_xticks([0, 17, len(gap.kpoints)1]) ax['gap0'].set_xticklabels(['0.5H', 'A', '0.5L']) plt.savefig('tz2_2_bandstructure.pdf') plt.show()
Run it using:
python ../tz2_2_bandstructure.py
You should get a less resolved version of the following figure:
The path of this band structure is first along a highsymmetry path in the k_z=0 plane (M\GammaK). Then is goes out of plane to the k_z=\pi/c plane (KH). It then goes in along the same path as in the k_z=0 plane, but translated to the k_z=\pi/c one (HAL).
At first glance, the two figures in the right column look very similar, as they are both gapped. The band gap energy is slightly smaller at 5 GPa compared to the 0 GPa value, which seems to indicate that the band gap of BiTeI closes with increasing pressure. However, if we compute the pressure dependence of the band gap energy between these two pressures, we will observe that the band gap decreases up to P\sim 2 GPa (and almost vanishes), then starts to \textit{increase} as pressure is further increased (see Figure 4 of [BrousseauCouture2020]). Knowing that a topological phase transition must be accompanied by a closing of the band gap, this behavior suggests that one of those pressures might have a nontrivial topology. If this is the case, the electronic bands in the nontrivial phase would display a band inversion.
Around the A point, one can notice some band dispersions characteristic of band inversion. In order to check whether or not there is inversion of character between bands which are gapped due to spinorbit coupling, we plot the fatbands. One can use the fatbands plotting feature of abipy as used in the Wannier90 tutorial, but instead here we will use a more advanced script which plots the difference in character between the p orbitals of bismuth and of those of both the tellurium and iodine. This procedure makes more clear the band inversion and the script is the following:
#!/usr/bin/env python from abipy.abilab import abiopen import abipy.data as abidata import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.axes_grid1 import make_axes_locatable def get_projection(fb, l, typat, blist=None, fact=1.0): # This function was adapted from the plot_fatbands_typeview function e0 = fb.ebands.get_e0('fermie') x = np.arange(fb.nkpt) mybands = range(fb.ebands.mband) if blist is None else blist eigens = np.zeros((fb.nkpt, len(mybands))) width = np.zeros_like(eigens) for symbol in typat: for spin in range(fb.nsppol): wl_sbk = fb.get_wl_symbol(symbol) * (fact / 2) for ib, band in enumerate(mybands): eigens[:, ib] = fb.ebands.eigens[spin, :, band]  e0 w = wl_sbk[l, spin, band] width[:, ib] += w return eigens, width def plot_relative_projection(eigen, proj1, proj2, ax, x): # comptute relative projection and transform to [0,1] interval diff = proj1  proj2 diff = np.amin(diff) maxproj = np.amax(diff) diff = diff/maxproj for ib in range(np.shape(diff)[1]): ax0 =ax.scatter(range(x), eigen[:, ib], c=diff[:, ib], marker='o', s=80, cmap='seismic', edgecolor='k', zorder=2, vmin=0, vmax=1) return ax0 flist = ["tz2_2o_DS13_FATBANDS.nc", "tz2_2o_DS23_FATBANDS.nc"] title = ['0GPa', '5GPa'] bands = [np.arange(36,40), np.arange(36,40)] fig, ax = plt.subplots(1, len(flist), figsize=(5*len(flist), 6)) for j in range(len(flist)): # plot electronic bands data = abiopen(flist[j]) data.ebands.set_fermie_to_vbm() data.ebands.plot(with_gaps=False, ylims=(0.4, 0.8), ax=ax[j], show=False, fontsize=16, zorder=1) # compute Bip projections bi_eigen, bi_proj = get_projection(data, l=1, typat=['Bi'], blist=bands[j]) # compute the sum of Tep and Ip projections tei_eigen, tei_proj = get_projection(data, l=1, typat=['Te', 'I'], blist=bands[j]) # compute and plot relative projection relplot = plot_relative_projection(bi_eigen, bi_proj, tei_proj, x=data.nkpt, ax=ax[j]) ax[j].set_title(title[j]) ax[j].set_xticks([0, 17, data.nkpt1]) ax[j].set_xticklabels(['0.5H', 'A', '0.5L'], fontsize=12) divider = make_axes_locatable(ax[1]) cax = divider.append_axes("right", size="7%", pad=0.15) cbar = plt.colorbar(relplot, cax=cax, ticks=[0, 1]) cbar.ax.set_yticklabels(['(Te+I)$_p$', 'Bi$_p$']) plt.savefig('tz2_2_relative_projection.png') plt.show()
cp r ../input_tz2_3/ .
python ../tz2_2_relative_projection.py
which produces a less converged version of this figure:
At 0 GPa, there is no band inversion while at 5 GPa, there is one around the A point. We thus suspect that there is a nontrivial Z2 invariant for this material, that we will now compute using Z2Pack.
3. Computing the Z2 invariant using Z2Pack¶
You need first to install Z2Pack, if not yet available on your computer. Please, see Z2Pack home, or directly the Installation and setup information.
Then, to run the example presented in this section, you need to download four big files from GitHub, https://github.com/abinit/abinit_assets/tree/main/tests/tutoplugs/Input/results_tz2_3 , and store them in the directory results_tz2_3, such that
ls results_tz2_3
BiTeI_0.msgpack BiTeI_1.msgpack BiTeI_2.msgpack BiTeI_3.msgpack
You will now launch the script for ABINIT with Z2Pack:
cp r ../input_tz2_3/ .
python ../z2.py
given by the following file:
#!/usr/bin/env python # * coding: utf8 * # Z2 invariant calculation for BiTeI at 0 and 5 GPa import os import z2pack if not os.path.exists('./results_tz2_3'): os.makedirs('./results_tz2_3') # creating the System object # We will create one for each pressure. BiTeI_0gpa = z2pack.fp.System( input_files=['input_tz2_3/tz2_3_0gpa.abi', 'input_tz2_3/wannier90.win' ], kpt_fct=z2pack.fp.kpoint.abinit, kpt_path='tz2_3_0gpa.abi', command='abinit tz2_3_0gpa.abi >& log', executable='/bin/bash' ) BiTeI_5gpa = z2pack.fp.System( input_files=['input_tz2_3/tz2_3_5gpa.abi', 'input_tz2_3/wannier90.win' ], kpt_fct=z2pack.fp.kpoint.abinit, kpt_path='tz2_3_5gpa.abi', command='abinit tz2_3_5gpa.abi >& log', executable='/bin/bash' ) # calculating the WCC # Invariant at fixed k_z planes result_0 = z2pack.surface.run(system=BiTeI_0gpa, surface=lambda s, t: [s / 2, t,0.0], num_lines=20, save_file = './results_tz2_3/BiTeI_0.msgpack', load=True) result_1 = z2pack.surface.run(system=BiTeI_0gpa, surface=lambda s, t: [s / 2, t,0.5], num_lines=20, save_file = './results_tz2_3/BiTeI_1.msgpack', load=True) result_2 = z2pack.surface.run(system=BiTeI_5gpa, surface=lambda s, t: [s / 2, t,0.0], num_lines=20, save_file = './results_tz2_3/BiTeI_2.msgpack', load=True) result_3 = z2pack.surface.run(system=BiTeI_5gpa, surface=lambda s, t: [s / 2, t,0.5], num_lines=20, save_file = './results_tz2_3/BiTeI_3.msgpack', load=True) print('for 0 GPa:') print('z2 topological invariant at kz = 0.0: {0}'.format(z2pack.invariant.z2(result_0))) print('z2 topological invariant at kz = 0.5: {0}'.format(z2pack.invariant.z2(result_1))) print('for 5 GPa:') print('z2 topological invariant at kz = 0.0: {0}'.format(z2pack.invariant.z2(result_2))) print('z2 topological invariant at kz = 0.5: {0}'.format(z2pack.invariant.z2(result_3)))
It works as follow:

We define the system using z2pack.fp.System().

We define the surfaces for which the invariant is computed using z2pack.surface.run. We need two timereversal invariant planes, in this case the k_z=0 and k_z=\pi/c planes. This function also runs the calculation. It will create a directory where for each new line on the surface, corresponding to different points in the x direction, a computation using ABINIT+Wannier90 is performed. This data is preserved and if you want to add a line, you can use the argument load=True which will restart the calculation and perform only the lines that were not already done. This is what is done in this tutorial, as in order for this calculation to be sufficiently converged, it takes too much time.

The surface invariant is stored in the z2pack.invariant.z2 variable.
The script delivers:
for 0 GPa:
z2 topological invariant at kz = 0.0: 0
z2 topological invariant at kz = 0.5: 0
for 5 GPa:
z2 topological invariant at kz = 0.0: 0
z2 topological invariant at kz = 0.5: 1
which informs us that at 0 GPa, the system is a trivial insulator, while at 5 GPa, it is topologically nontrivial.
We included a script that helps visualize the trajectories of the Wannier charge centers. Run it like this:
python ../tz2_3_hwcc.py
You will find the following figure:
In this figure each open circle represents a hybrid Wannier charge center (HWCC). The left column is calculated in the k_z=0 plane and the right column in the k_z=0.5 plane. The HWCC form “bands”, similar to the band structure of a dispersion relation, but the values are defined mod 1, so the y axis is periodic. In the continuum case, the parity of the number of times a horizontal line crosses the HWCCs constitutes a topological invariant of the k_z=0 and 0.5 planes, which we denote by \Delta. If \Delta differs between k_z=0 and k_z=0.5, the bulk Z_2 invariant is \nu=1.
To calculate this invariant on a coarse kmesh, we consider the points marking the middle of the largest gap at each k_x value rather than a horizontal line. These points are marked by blue diamonds above. Whenever the location of the middle of the gap changes between two adjacent k_x values, we count the number of HWCC that exist between the two gap centers and sum this number for all the crossings as k_x goes from 0 to \pi/a_x. If this value is even, \Delta=0, and if it is odd, \Delta=1. Performing this counting procedure for both k_z planes produces \Delta=0 for each case besides k_z=0.5 with P=5 GPa, indicating that the 0 Gpa phase is trivial and the 5 GPa phase is topological.
4. Conclusion¶
We have seen that in material with strong spinorbit coupling, band crossing at the Fermi level can lead to band inversion seen in the orbital character of the bands. The gap thus created can be topologically nontrivial. To know whether or not a material is a topological insulator in the sense of the Z2 invariant, we can look at the displacement of the Wannier charge centers instead of integrating the whole Berry curvature. This leads to the calculation performed by Z2Pack.
For BiTeI, we found that at 0 GPa, the material is a nontrivial insulator. Applying pressure, the bands reach a point where they cross at the Fermi level. Upon increasing pressure, the bands countinue to move and spinorbit ocupling opens a nontrivial topological gap, which makes this material a topological insulator. With the help of firstprinciples calculation using ABINIT, along with clever ways to compute the topological invariant using Z2Pack, we were able to confirm this intrinsic property of BiTeI.