3D topological insulator¶
Note
You need to use the newer xfbp/xfplo version,
which comes with pyfplo
in order for the cmd
scripts to work
properly.
Here we discuss the calculation of the topological indices for 3d topological insulators. The algorithm is based on Wannier centers as in 2D topological insulator. The only difference is that we need to calculate four \(Z_2\) invariants for 4 planes in the 3d BZ. In principle there are 6 different planes forming the faces of a parallelepiped whose corners are the eight distinct time reversal invariant momenta (TRIM) of any chosen smallest primitive reciprocal cell. However, if the electronic spectrum is gapped above a certain band there are only 4 independent \(Z_2\) invariants, which are usually grouped into four topological indices \(\nu_0;(\nu_1\nu_2\nu_3)\). \(\nu_0=1\) indicates a strong topological insulator while any \(\nu_{1,2,3}=1\) indicate conditions on the number of surface state pairs for certain surface orientations. If \(\nu_0=0\) but some \(\nu_{1,2,3}\neq 0\) it is a weak topological insulator. If all indices are zero it is a trivial insulator. The mapping of \(Z_2\) invariants to indices is as follows: if any two parallel planes of the parallelepiped have differing invariants \(\nu_0=1\). We take the two planes spanned by the first two reciprocal lattice vectors (\(g_{1,2}\)) called \(z_0\) if the plane goes through the origin and \(z_1\) if it goes through \(g_{3}\). The three weak indices are identical to the \(Z_2\) invariants of the three planes, which do not go through the origin, which we call \(x_1\) (spanned by \(g_{2,3}\) through \(g_1\)), \(y_1\) (spanned by \(g_{3,1}\) through \(g_2\)) and \(z_1\) (spanned by \(g_{1,2}\) through \(g_3\)). The weak indices of course depend on the chosen TRIM points.
The algorithm selects the TRIM points and performs four \(Z_2\)
calculations. It tries to determine the invariants via the automatic
procedure explained in 2D topological insulator. The results including
reliability and electronic gap estimates are printed to the output and
convenience files Z2_3dTI_homo...cmd
for use in
xfbp are created.
This algorithm works for centro symmetric and non-centro symmetric
compounds. A version of this algorithm is linked directly into
fplo. The difference is that when used from fplo a all-orbital
Wannier basis is created instead of the reduces Wannier models, which
usually are the basis of pyfplo.slabify
. This leads to differently
locking Wannier center curves, since there are more orbitals (semi
core and such). If surface state calculations are planned it is anyway
necessary to create a Wannier model first. The resulting
\(Z_2\)-calculations are faster than in fplo.
For illustration we chose a system, which is centro symmetric, since
in this case we know the invariants from the parity algorithm, which
allows comparison. This example also shows some fo the complications.
We will use both the fplo version and the pyfplo
version of the
Wannier center algorithm to understand the pecularities.
The tutorial files are in
FPLO.../DOC/pyfplo/Examples/slabify/3dTI/Bi2Se3
.
To start with change into this directory and have a look into the topological insulator submenu in fedit. You will see that “Force wannier center” is switched on, to force the use of this algorithm despite the fact that this compound has an inversion center.
Now, have a look at makewandef.py
. It illustrates how to
create a =.wandef
quickly. There is a particularity in the
system, which is that the Bi6p and Se4p orbitals are separated from
lower and higher bands through two gaps. This makes creating a Wannier
model especially easy. We use the option lbands and ubands to
define an energy window which excludes all other bands.
Please execute
python makewandef.py
You should have a =.wandef file now. Next run fplo. E.g.:
fplo.... > out
to produce a bandstructure, +wancoeff
and the output from the
fplo Wannier center calculation. This will take some time, since the
fplo version of the WC algorithm is not the fastest, due to the
larger number of bands needed. Output related to the TI calculation
can be found via:
grep TI: out
First, the invariants from parities are printed and then the Wannier center algorithm is executed. Various files are created. Note, that the indices from parities can differ from the Wannier center results; first, because the WC algorithm is tricky or second, if the electronic gap for a certain homo is zero, in which case both algorithms makes no sense.
You will see the following output from the Wannier centers:
TI: Invariants:
TI: homo invariants E(k=0) estim.gaps [eV] for
TI: z0 z1 x1 y1
TI: 114 ? 0;(000) -1.09567 0.05585 0.01917 0.01917 0.01917
TI: 116 ? 1;(111) -0.79835 0.01306 0.03873 0.03873 0.03873
TI: 118 * 1;(000) -0.22789 0.44699 0.43513 0.43513 0.43513
TI: 120 ? 1;(111) 0.28647 0.15850 0.14845 0.14845 0.14845
TI: 122 ? 0;(111) 1.42635 0.02117 0.17378 0.17378 0.17378
TI: 124 ? 0;(000) 1.44752 0.21087 0.29692 0.29692 0.29692
TI:
TI: homo invariants E(k=0) Z2 (reliability)
TI: z0 z1 x1 y1
TI: 114 ? 0;(000) -1.09567 0 ( 77%) 0 ( 77%) 0 ( 77%) 0 ( 77%)
TI: 116 ? 1;(111) -0.79835 0 ( 50%) 1 ( 73%) 1 ( 73%) 1 ( 73%)
TI: 118 * 1;(000) -0.22789 1 (100%) 0 (100%) 0 (100%) 0 (100%)
TI: 120 ? 1;(111) 0.28647 0 ( 92%) 1 ( 65%) 1 ( 65%) 1 ( 65%)
TI: 122 ? 0;(111) 1.42635 1 ( 69%) 1 (100%) 1 (100%) 1 (100%)
TI: 124 ? 0;(000) 1.44752 0 ( 80%) 0 (100%) 0 (100%) 0 (100%)
The first table shows homos, indices the energy at the Gamma point and estimates of the electronic gap (taken from the values on the finite grid used for the WC algorithm), while the second shows the \(Z_2\) invariants and their estimated reliability for the four planes of the parallelepiped mentioned above. If the reliability is less than 100% a question mark is printed after the homo. The homo, which is likely the highest occupied band is marked by a * for orientation. This is of course only reliable if the system has a true gap.
In our case we see that there are rather small estimated gaps for lots of bands. Fortunately above the Fermi energy (homo 118) the gap is sizable. Not, surprisingly the reliabilities are all very high (actually 100%). We also find lots of question marks. Let’s have a look at the fplo and Wannier model bands (red/green).
We see the two gaps, separating the p-bands as discussed above and that nearly all bands are separated by gaps (although small). This means that topological invariants have a meaning. Note, that it is not needed that the gap is visible in the DOS. All one needs is that the bands of homo and homo+1 never cross (warped gap).
Since this table is not ment for straight forward consumption we have to inspect the Wannier centers by hand. Run:
xfbp Z2_3dTI_homo118.cmd
and so forth. You will get several pictures. Let’s start with Wannier centers and reference line for homo 118. As already pointed out in the table his band has a nice electronic gap above itself, happens to be the last occupied band below the Fermi level and has high reliability of the automatically determined \(Z_2\) invariants (reference line winding number).
It can be clearly seen that the \(x_1\), \(y_1\) and \(z_1\) planes have zero \(Z_2\) invariant, while plane \(z_0\) is non-trivial. To really verify this use the right mouse click in xfbp close to the higest Wannier center curve in the \(z_0\) panel close to \(\theta=1\) and \(k_y=0\). It will show that only one center (called Set118 in xfbp) comes down from it’s degenerate value of +-1. A clean straight reference line can be drawn for \(\theta=0.9\) (not shown), which only crosses this center and hence crosses an odd number of Wannier centers, which results in \(Z_2=1\). This and the fact that we can visually connect the Wannier center curves in a reasonable smooth way convinces us that the topological indices are 1;(000). (Remember \(z_0\neq z_1 \to \nu_0=1\) and \(\nu_{1}=Z_2(x_1), \ldots\).) In fact due to symmetry the tree planes \(x_1\), \(y_1\) and \(z_1\) give the same result in this compound.
Next we have a look at all the other homos, just for learning purposes.
The Wannier centers and reference line for homo 114 are quite messy. However, you should convince yourselfs that \(z_1\) has a region (\(\theta\approx0.5\)) without any centers crossing \(\to\) trivial, while for the \(z_0\)-plane an even number of curves cross any reference line. Hence we get 0;(000).
The Wannier centers and reference line for homo 116 are also messy and reveal another possible complication. First, note that both 114 and 116 yield curves with fast varying (steep slopes) curves. These steep curves generate only very few points on their steep section. Depending on the density of the ky mesh there might even be zero points on such sections, which would make it look as if certain curves end for some ky value to reemerge for a later ky-value. Any algorithm including our own judgment can fail in such cases. All we can do is to take a denser ky-mesh. In this case the automatic reference line algorithm is correct as can be visually verified. There is an even/odd number of curves crossed by a suitably chosen reference curve for the \(z_0\)/\(z_1\) planes (\(x_1,y_1=z_1\)). Hence we get 1;(111). Note, however, that the electronic gaps above homos 114 and 116 are rather small, they could well be zero. If this homo were important we would try to use much finer meshes to make absolutely sure that the gap is finite.
Homo 122 also has a small electronic gap. If we assume that it is finite we can proceed analysing the invariants. The Wannier centers and reference line for homo 120 and the Wannier centers and reference line for homo 122 both show that the indices printed in the output table are wrong (as can be verified in this case due to the parity algorithm).
For homo 120 the \(z_0\) plane is clearly trivial since we can draw a non crossing reference line at \(\theta=\pm1\) (the automatic blue reference line is just that). For the other three planes the algorithm made a mistake around ky=0.33, where a single red dot appears. This mistake is due to a steep section at one curve. One can find hand drawn reference curves, which cross an even number of times and hence these planes are trival. The result should be 0;(000). Please also note that the reliability for homo 120 is high for the \(z_0\) plane but only 65% for the other three, which fits to the findings explained above. This is of course only an indicator.
For homo 122 the situation is even trickier. According to the reliabilities the \(z_0\) panel might have a problem. Indeed, the \(z_1\) plane (and \(x_1\), \(y_1\)) is clearly topological, but the \(z_0\) plane was determined as non trivial although it is not. Again due to small electronic gaps in the \(z_0\) plane we get very steep curves close to ky=0. In fact this curve is so steep that we see only one point in the steep section. In xfbp use the right mouse click at the dot in the upper left corner. It will show that there are two curves. Given the rest of the dots it is clear that one dot is connected with the declining curve and the other must fall very steeply to connect to a curve close to \(\theta=0.28\). Hence the reference line actually crosses two curves and not one. Consequently, \(z_0\neq z_1 \to\) 1;(111) in contrast with the result printed in the output table.
Homo 124 is again a simple case and will not be discussed here.
In summary we saw that the fplo based version of the procedure is
rather slow, which is bad if a high grid density is needed to obtain
reliable results. Let’s continue with the pyfplo
version then. We
already set up the Wannier basis production in our first fplo
run. Now, we can proceed by re-running:
fplo...
Note that we do not use an output file (we don’t want to overright the
TI results). This will produce among others the file +hamdata
,
which we will use for the pyfplo
based version of the TI procedure.
If the fplo run is over we confirm the Wannier basis by checking:
xfbp wband.cmd
which should show green bands (finite cut-off tight-binding Wannier
bands) on top of red bands (exact Wannier bands) on top of black bands
(fplo bands). Please, change into slabify/Z2
and execute:
python -u Z2.py | tee out
The flag -u
un-buffers the output, such that the progress is updated
in real time. The file’s content is
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 | #! /usr/bin/env python
from __future__ import print_function
import sys,copy
import numpy as np
import pyfplo.slabify as sla
import pyfplo.common as com
import pyfplo.fploio as fploio
# ===================================================================
#
# ===================================================================
def work():
print( '\nversion: ',sla.version,'\n')
np.set_printoptions(precision=5,suppress=True,linewidth=120)
hamdata='../../+hamdata'
s=sla.Slabify()
s.object='3d'
s.printStructureSettings()
s.prepare(hamdata)
p=fploio.INParser()
p.parseFile('../../=.in')
d=p()('special_sympoints')
l=[]
for i in range(d.size()):
l.append([ d[i]('label').S,d[i]('kpoint').listD])
bp=com.BandPlot()
bp.points=l
bp.ndiv=100
bp.calculateBandPlotMesh(s.dirname)
#help(sla)
s.calculateBandStructure(bp)
# 14 -> 0;(000)
# 16 -> 1;(111)
# 18 -> 1;(000)
# 20 -> 0;(000)
# 22 -> 1;(111)
# 24 -> 0;(000)
s.calculate3dTIInvariants(200,200,homos=range(14,25,2),efhomo=18)
return
# ===================================================================
#
# ===================================================================
if __name__ == '__main__':
work()
|
In lines 34-47 we calculate the band structure of the Wannier model along the same high symmetry points as were used in the fplo calculation. With the help of:
xfbp slabifyres/+band_sf
and a right mouse click we find out that the Fermi level is above homo
18 (just a coincidence that 118-18=100). This is used in line 56 to
define a set of homos and efhomo. The grid is chosen finer for
both the “integration” and the ky direction. See calculate3dTIInvariants
.
The output table:
Invariants:
homo invariants E(k=0) estim.gaps [eV] for
z0 z1 x1 y1
14 ? 0;(000) -1.09645 0.02338 0.01563 0.01563 0.01563
16 ? 1;(111) -0.79820 0.00531 0.02098 0.02098 0.02098 <-small
18 * 1;(000) -0.22063 0.42925 0.43611 0.43611 0.43611
20 0;(000) 0.27854 0.14357 0.06937 0.06937 0.06937
22 1;(111) 1.42038 0.00911 0.17176 0.17176 0.17176 <-small
24 0;(000) 1.42948 0.20549 0.24304 0.24304 0.24304
homo invariants E(k=0) Z2 (reliability)
z0 z1 x1 y1
14 ? 0;(000) -1.09645 0 (100%) 0 ( 93%) 0 ( 93%) 0 ( 93%)
16 ? 1;(111) -0.79820 0 (100%) 1 ( 86%) 1 ( 93%) 1 ( 93%)
18 * 1;(000) -0.22063 1 (100%) 0 (100%) 0 (100%) 0 (100%)
20 0;(000) 0.27854 0 (100%) 0 (100%) 0 (100%) 0 (100%)
22 1;(111) 1.42038 0 (100%) 1 (100%) 1 (100%) 1 (100%)
24 0;(000) 1.42948 0 (100%) 0 (100%) 0 (100%) 0 (100%)
looks a bit different from the previous one. First of all the gap values are different since there are always tiny shifts in a reduced Wannier model compared to a full band structure and second we used a finer grid. Furthermore, the reliabilities increased considerably, again mostly due to the finer grid. Finally, all topological indices are correct.
We can now load all the Z2_3dTI_homo...cmd
files into xfbp
and analyse the validity by hand. We will look at three Wannier
center graphs. The Wannier centers and reference line for homo 16 show a quick change
in the \(z_0\)-panel around ky=0.16, which is most likely due to the small
electronic gap. Furthermore, the steep sections in the \(z_1\)-panel are
now represented by more points, which makes the automatic reference
line algorithm more reliable.
The Wannier centers and reference line for homo 20 has the correct reference line due to the increased mesh size.
The Wannier centers and reference line for homo 22 now clearly shows that two WC curves meet close to the upper left corner of the \(z_0\)-panel, which makes this plane trivial. Consequently, the invariants are correctly determined as 1;(111).
In summary, the use of a reduce Wannier function model leads to a much faster TI procedure, which allows to use finer grids, which in turn increase the accuracy of the results. Additionally, one could calculate surface spectra using the same model.