2.3. Capabilities of WannierTools

• Bulk band calculation (points mode, line mode and plane mode)
• BulkFS calculation
• BulkFS plane calculation
• Bulk spin texture calculations
• Density state(DOS) calculations
• Find Nodes calculation
• Energy gap calculations (plane and cube mode)
• Slab band calculation
• Nanowire/nanoribbon band calculation
• Surface state ARPES calculation
• Surface state QPI calculation
• Fermi arc calculation
• Spin texture calculation
• Berry phase calculation
• Berry curvature calculation for 3D bulk case
• Berry curvature calculation for slab system
• Anomalous Hall conductivity (AHC)
• Wannier charge center/Wilson loop calculation
• Z2 number for 3D bulk materials
• Chern number for 3D bulk materials
• Mirror Chern number calculation New
• Weyl Chirality calculation
• Landau level calculations
• Ordinary magnetoresistance calculations New

2.3.1. Bulk band calculation (points mode, line mode and plane mode)

Points mode

You can calculate the energy bands with the given k points in the KPOINTS_3D KPOINTS_3D card.

Input

Typical flags for this mode in the wt.in.

&CONTROL
BulkBand_points_calc = T
/

KPOINTS_3D
4              ! number of k points
Direct         ! Direct or Cartesian
0.00000 0.00000 0.0000
0.00000 0.00000 0.5000
0.50000 0.50000 0.0000
0.00000 0.00000 0.0000

Output

The outputs for this mode is bulkek-pointsmode.dat

Line mode

Calculate bulk energy band for a series k lines. This is the basic calculation after the construction of Wannier functions. You have to compare your Wannier interpolated bands with the DFT bands. Those two bands should match well around the Fermi level.

Input

Typical flags for bulk band calculation in the wt.in.

&CONTROL
BulkBand_calc = T
/
&PARAMETERS
Nk1 = 101   ! Number of k points for each k line
/

KPATH_BULK     ! k point path
4              ! number of k lines only for bulk band
G 0.00000 0.00000 0.0000 Z 0.00000 0.00000 0.5000
Z 0.00000 0.00000 0.5000 F 0.50000 0.50000 0.0000
F 0.50000 0.50000 0.0000 G 0.00000 0.00000 0.0000
G 0.00000 0.00000 0.0000 L 0.50000 0.00000 0.0000

See CONTROL, PARAMETERS, KPATH_BULK

Output

The outputs for bulk band calculation are bulkek.dat and bulkek.gnu. You can get the band plot by running

gnuplot bulkek.gnu

or

xmgrace bulkek.dat

to get bandstucture plot.

The data structure for bulkek.dat

0.000000000       -2.673821992  119   80   80  119   80   80  205  138  138   70   40   40   70   40   40    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
0.016453872       -2.681536808  118   78   78  118   78   78  203  134  134   82   41   41   82   41   41    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0

1. The 1st column represents k points for the given kpath (KPATH_BULK)
2. The 2nd column is the energy level
3. From the 3rd to the n’th column are the projected weight of the wave function at each k point and each band onto each wannier orbitals. Those weights are normalized to 255 for the color plot convinence.

The subrotine for this feature is ek_bulk.f90 .

Plane mode

Calculate band structure in a k slice(plane) specified by KPLANE_BULK card. The mode is very useful to visualize the Dirac/Weyl cone. You have to set the following tags in wt.in

&CONTROL
BulkBand_plane_calc = T
/
&PARAMETERS
Nk1 = 51   ! Number of k points along the first vector in KPLANE_BULK
Nk2 = 51   ! Number of k points along the second vector in KPLANE_BULK
/

KPLANE_BULK   ! fractional coordinates
 0.00  0.00  0.30   ! Middle point for a k slice(plane) in 3D BZ. Usually, the position of Dirac points.
 0.50  0.00  0.00   ! The first vector to define k plane(slice) in 3D BZ
 0.00  0.50  0.00   ! The second vector to define k plane(slice) in 3D BZ

The output file is bulkek_plane.dat, bulkek_plane-matlab.dat and bulkek_plane.gnu. You can get bulkek_plane.png with

gnuplot bulkek_plane.gnu

The bulkek_plane-matlab.dat is in MATLAB data format. You can plot the Dirac cone with matlab.

The format of bulkek_plane.dat is as follows:

        # kx                 ky                 kz                 k1                 k2                 k3   E(Numoccupied-1)     E(Numoccupied)   E(Numoccupied+1)   E(Numoccupied+2)
-0.299354337       -0.518496963        0.180167841       -0.518496936       -0.299354384        0.180167841       -1.311721381       -1.311705191        0.588683811        0.588872215
-0.299354337       -0.504670376        0.180167841       -0.511583643       -0.287380208        0.180167841       -1.294078082       -1.293904952        0.586780093        0.587249790
...

Column 1-3rd are k points in cartesian coordinates. Column 4-6th are k points in a rotated cartesian coordinates where the x and y direction are line in the k plane and the z direction is perpendicular to the k plane you specified. Column 7-10th are energies at each k point. Here we only print out 4 energy bands around the fermilevel. It depends on NumOccupied. Usually, I choose column 4th and 5th as k coordinates and choose 8 and 9 as energy bands to show the Dirac cone shown below.

2.3.2. BulkFS calculation

Bulk Fermi surface calculation.

Input

You should specify the number of k points for each three reciprocal vectors Nk1, Nk2, Nk3 in NAMELISTS PARAMETERS

&CONTROL
BulkFS_calc = T
/
&PARAMETERS
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
Nk2 = 101   ! No. of slices for the 2nd reciprocal vector
Nk3 = 101   ! No. of slices for the 3rd reciprocal vector
/

KCUBE_BULK
  0.00  0.00  0.00   ! Original point for 3D k plane
  1.00  0.00  0.00   ! The first vector to define 3d k space plane
  0.00  1.00  0.00   ! The second vector to define 3d k space plane
  0.00  0.00  1.00   ! The third vector to define 3d k cube

See CONTROL, PARAMETERS

Output

The outputs for this function are FS3D.bxsf. You can plot the FS with xcrysden run

xcrysden --bxsf FS3D.bxsf

to get the plot.

By the way, Bulk band and BulkFS calculations were already implemented in Wannier90 code.

2.3.3. BulkFS plane calculation

Bulk Fermi surface in a fixed k plane specified by KPLANE_BULK

Input

You should specify the number of k points for each three reciprocal vectors Nk1, Nk2 in NAMELISTS PARAMETERS

&CONTROL
BulkFS_Plane_calc = T
/
&PARAMETERS
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
Nk2 = 101   ! No. of slices for the 2nd reciprocal vector
/

KPLANE_BULK  ! in fractional coordinates
  0.00  0.00  0.00   ! Original point for 3D k plane
  1.00  0.00  0.00   ! The first vector to define 3d k space plane
  0.00  1.00  0.00   ! The second vector to define 3d k space plane

See CONTROL, PARAMETERS

Output

The outputs for this function are fs.gnu, fs.png.

gnuplot fs.gnu

to get the plot.

2.3.4. Bulk spin texture calculations

Calculate spin texture for bulk system that with vacuum or without inversion symmetry. For the bulk system with vacuum, you can calculate the surface projected spin texture. This is useful for comparing with the ARPES experiments. if you cut a slab system from a periodic tight binding model, then there is no charge relaxation on the surface which would change the surface state a lot. In this case, you have to do the first-principle calculations for a finite thickness slab system that with vacuum. Then you can construct Wannier functions for this system and use WannierTools to get the iso-energy plot of the Fermi surface (BulkFS_plane_calc =T) and get the surface projected spin texture (Bulkspintext_calc=T).

There is one example in the examples/Bi2Se3-6Qlayers.

2.3.5. Density state(DOS) calculations

Calculation density of state for the bulk system. The typical setup in wt.in:

&CONTROL
DOS_calc = T
/
&PARAMETERS
OmegaNum = 601    ! number of slices of energy
OmegaMin = -1.0   ! erergy range for DOS
OmegaMax =  1.0
Nk1 = 51   ! No. of slices for the 1st reciprocal vector
Nk2 = 51   ! No. of slices for the 2nd reciprocal vector
Nk3 = 51   ! No. of slices for the 3nd reciprocal vector
/

KCUBE_BULK
  0.00  0.00  0.00   ! Original point for 3D k plane
  1.00  0.00  0.00   ! The first vector to define 3d k space plane
  0.00  1.00  0.00   ! The second vector to define 3d k space plane
  0.00  0.00  1.00   ! The third vector to define 3d k cube

Outputs are dos.dat and dos.gnu. dos.eps will be obtained with

gnuplot dos.gnu

2.3.6. Energy gap calculations (plane and cube mode)

We support two modes for energy gap calculations.The formula is \(gap(k)= E_{NumOccupied+1}(k)- E_{NumOccpuied}(k)\)

Gap_plane mode

Calculate the energy gap for the k points in the KPLANE_BULK. This is useful to show Weyl points and nodal line structure.

Input

Typical input parameters for BulkGap_plane calculation

&CONTROL
BulkGap_Plane_calc = T
/
&PARAMETERS
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
Nk2 = 101   ! No. of slices for the 2nd reciprocal vector
/

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

See CONTROL, PARAMETERS, KPATH_BULK

Output

The outputs for Gap_plane mode are GapPlane.dat, GapPlane.gnu. The GapPlane.png will get by

gnuplot GapPlane.gnu

The head of GapPlane.dat

kx              ky              kz             gap             Ev4             Ev3             Ev2             Ev1             Ec1             Ec2             Ec3             Ec4              k1              k2              k3
0.00000000      0.00000000      0.00000000      0.45569845     -0.69109275     -0.69109055     -0.29654328     -0.29654073      0.15915772      0.15915871      1.24348171      1.24348457      0.00000000      0.00000000      0.00000000
0.03796028     -0.02191637      0.00548462      0.43770730     -0.77636510     -0.77598312     -0.26035113     -0.26027881      0.17742849      0.17771545      1.29499437      1.29505298      0.00000000      0.02500000      0.00000000

• Column 1-3 are the Cartesian coordinates of the k points in the KPLANE_BULK, in unit of \(\frac{1}{Angstrom}\)
• Column 4 is the energy gap
• Column 5-12 are the energy value for valence and conduction bands close to the Fermi level. There are 4 conduction bands and 4 valence bands.
• Column 13-15 are the Direct coordinates of the k points in the KPLANE_BULK

Gap_Cube mode

This helps us to find Weyl points and nodal line structure in the 3D BZ.

Input

Typical input parameters for BulkGap_cube calculation

&CONTROL
BulkGap_Cube_calc = T
/
&PARAMETERS
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
Nk2 = 101   ! No. of slices for the 2nd reciprocal vector
Nk3 = 101   ! No. of slices for the 3rd reciprocal vector
/

KCUBE_BULK
-0.50 -0.50 -0.50   ! Original point for 3D k plane
 1.00  0.00  0.00   ! The 1st vector to define 3d k cube
 0.00  1.00  0.00   ! The 2nd vector to define 3d k cube
 0.00  0.00  1.00   ! The 3rd vector to define 3d k cube

See CONTROL, PARAMETERS, KCUBE_BULK

Output

The outputs for Gap_plane mode are GapCube.dat, GapCube.gnu. The GapCube.png will get by

gnuplot GapCube.gnu

The head of GapCube.dat are

kx (1/A)        ky (1/A)        kz (1/A)      Energy gap              Ev              Ec      k1 (2pi/a)      k2 (2pi/b)      k3 (2pi/c)
0.00000000      0.87665487     -0.54846229      0.79075142     -0.34827281      0.44247861     -0.50000000     -0.50000000     -0.50000000
0.00000000      0.87665487     -0.51555455      0.86792416     -0.38635069      0.48157346     -0.50000000     -0.50000000     -0.45000000

• Column 1-3 are the Cartesian coordinates of the k points where energy gap is small than Gap_threshold, in unit of \(\frac{1}{Angstrom}\)
• Column 4 is the energy gap. Those values are smaller than Gap_threshold, see PARAMETERS
• Column 5-6 are the energy value for valence and conduction bands close to the Fermi level. There are 4 conduction bands and 4 valence bands.
• Column 7-9 are the Direct coordinates of the k points.

2.3.7. Find Nodes calculation

Beside by using GapCube and GapPlane to find Weyl/Dirac nodes or node lines, we can directly using FindNodes function. \(gap(k)= E_{NumOccupied+1}(k)- E_{NumOccpuied}(k)\)

Input

Typical input parameters for FindNodes_cube calculation

&CONTROL
FindNodes_calc = T
/
&PARAMETERS
Nk1 = 8   ! No. of slices for the 1st reciprocal vector
Nk2 = 8   ! No. of slices for the 2nd reciprocal vector
Nk3 = 8   ! No. of slices for the 3rd reciprocal vector
Gap_threshold = 0.0001 ! a value to determine which point should be identified as a node
/

KCUBE_BULK
-0.50 -0.50 -0.50   ! Original point for 3D k plane
 1.00  0.00  0.00   ! The 1st vector to define 3d k cube
 0.00  1.00  0.00   ! The 2nd vector to define 3d k cube
 0.00  0.00  1.00   ! The 3rd vector to define 3d k cube

Note

Please don’t set Nk1, Nk2, Nk3 too large. Otherwise, it will become very time consuming. Usually, 15*15*15 is enough to get converged number of Weyl/Dirac points.

Output

Outputs are Nodes.dat and Nodes.gnu. Nodes.png will be obtained by

gnuplot Nodes.gnu

Here are heads of output for WTe2 Nodes.dat

# local minimal position and the related energy gap
#      kx          ky          kz         gap           E          k1          k2          k3
    0.219436   -0.045611   -0.000001    0.000000    0.056688    0.121432   -0.045363   -0.000003
   -0.219515   -0.045063   -0.000001    0.000000    0.056461   -0.121476   -0.044818   -0.000002
    0.220195   -0.038682   -0.000002    0.000000    0.051264    0.121852   -0.038472   -0.000003
   -0.220183   -0.038936   -0.000001    0.000000    0.051618   -0.121845   -0.038724   -0.000003
    0.219514    0.045063    0.000001    0.000000    0.056459    0.121475    0.044818    0.000003
   -0.219434    0.045620    0.000002    0.000000    0.056692   -0.121431    0.045371    0.000004
   -0.220194    0.038678    0.000000    0.000000    0.051259   -0.121851    0.038468    0.000001
    0.220181    0.038941    0.000000    0.000000    0.051620    0.121844    0.038729    0.000001

You will find that there are 8 Weyl points in the BZ as expected.

2.3.8. Weyl Chirality calculation

After you identify the positions of Weyl points, you could use this function to calculate the chirality, which tells you whether a Weyl point is a sink or a source of the Berry Curvature.

Input

Typical input parameters for WeylChirality_calc calculation

&CONTROL
WeylChirality_calc = T
/
&PARAMETERS
Nk1 = 41   ! No. of slices for the 1st reciprocal vector, berry phase integration direction
Nk2 = 21   ! No. of slices for the 2nd reciprocal vector
/

WEYL_CHIRALITY
8            ! Num_Weyls
Cartesian    ! Direct or Cartesian coordinate
0.004        ! Radius of the ball surround a Weyl point
 0.219436   -0.045611   -0.000000    ! Positions of Weyl points, No. of lines should larger than Num_weyls
-0.219515   -0.045063   -0.000000
 0.220195   -0.038682   -0.000000
-0.220183   -0.038936   -0.000000
 0.219514    0.045063    0.000000
-0.219434    0.045620    0.000000
-0.220194    0.038678    0.000000
 0.220181    0.038941    0.000000

Output

Outputs are wanniercenter3D_Weyl.dat and wanniercenter3D_Weyl_i.gnu. wanniercenter3D_Weyl.png will be obtained by

gnuplot wanniercenter3D_Weyl_i.gnu

for ((i=1; i<9; i++)); do gnuplot wanniercenter3D_Weyl_$i.gnu;done

Note

i is an integer from 1 to Num_weyls

Here are heads of output for WTe2 wanniercenter3D_Weyl.dat

# Chirality              -1               1               1              -1               1              -1               1              -1
       # k            phase           phase           phase           phase           phase           phase           phase           phase
 0.00000000      0.99970932      0.00005854      0.00004671      0.99975139      0.00005851      0.99970861      0.00004736      0.99975087
 0.05000000      0.89229069      0.08696587      0.08941971      0.90855415      0.08723118      0.89170870      0.09022452      0.90795187
 0.10000000      0.79659821      0.16589558      0.17112299      0.82248889      0.16697194      0.79511289      0.17279423      0.82108022

The first line shows the chiralities of each Weyl point. The first column is k point. From the 2nd to the last column show the Wannier charge center phase. In total, there are Num_weyls columns.

2.3.9. Slab band calculation

Before using iterative green’s function to get the surface state spectrum for semi-infinite system. We also can just construct a finite thickness slab system and calculate the band structure for it.

Note

For slab calculations, please read carefully the input card SURFACE

Input

&CONTROL
SlabBand_calc = T
/
&SYSTEM
NSLAB = 10
/
&PARAMETERS
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
/
KPATH_SLAB
2        ! numker of k line for 2D case
K 0.33 0.67 G 0.0 0.0  ! k path for 2D case
G 0.0 0.0 M 0.5 0.5

See CONTROL, SYSTEM PARAMETERS, KPATH_SLAB

Output

Outputs are slabek.dat and slabek.gnu

The heads of slabek.dat are

0.0000000     -4.9575466     240
0.0508687     -5.0110528     226
0.1017373     -5.0566963     221
0.1526060     -5.0671994     220
...

• The 1st column are k points in the KPATH_SLAB
• The 2nd column are energy values.
• The 3rd column represent the surface weight, which is normalized to 255.

The colorfull plot slabek.png of the slab energy bands can be obtained by

gnuplot slabek.gnu

2.3.10. Nanowire/nanoribbon band calculation

Band calculation for wire system. Only one direction is periodic, the other two directions are confined.

Input

You don’t have to set the k path, because it only has one direction.

&CONTROL
WireBand_calc = T
/
&SYSTEM
NSLAB1 = 4
NSLAB2 = 4
/
&PARAMETERS
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
/

See CONTROL, SYSTEM PARAMETERS

Output

Outputs are ribbonek.dat and ribbonek.gnu. The data format of ribbonek.dat is the same as slabek.dat. Get plot ribbonek.png with

gnuplot ribbonek.gnu

2.3.11. Surface state ARPES calculation

One important feature for topological materials is the surface state. The bulk-edge correspondence tells us, if the topological property of the bulk system is nontrivial, then there will be nontrivial states on the surface. Nowadays, there are several method to detect the surface states. One direct way is the angle resolved photoemission spectroscopy (ARPES). Such spectrum can be obtained by the iterative Green’s function.

Note

For slab calculations, please read carefully the input card SURFACE

Input

&CONTROL
SlabSS_calc = T
/
&PARAMETERS
OmegaNum = 101
OmegaMin = -1.0
OmegaMax =  1.0
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
NP = 2      ! principle layer
/
KPATH_SLAB
2        ! numker of k line for 2D case
K 0.33 0.67 G 0.0 0.0  ! k path for 2D case
G 0.0 0.0 M 0.5 0.5

See CONTROL, PARAMETERS, KPATH_SLAB

NP : integer valued, Number of principle layers, see details related to iterative green’s function. Used if SlabSS_calc= T, SlabArc_calc=T, SlabSpintexture_calc=T. Default value is 2. You need to do a convergence test by setting Np= 1, Np=2, Np=3, and check the surface state spectrum. Basically, the value of Np depends on the spread of Wannier functions you constructed. One thing should be mentioned is that the computational time grows cubically of Np.

Output

The output files are surfdos_l.dat, surfdos_r.dat, surfdos_l.gnu, surfdos_r.gnu. You could get the the spectral function plots with

gnuplot surfdos_l.gnu
gnuplot surfdos_r.gnu

_l and _r means the top and bottom surface.

2.3.12. Surface state QPI calculation

Settings for this feature are almost the same as Fermi arc calculation. Only difference is that you should set

# please set SlabQPI_kplane_calc = T from V2.4.2
&CONTROL
SlabQPI_kplane_calc          = T
/

Output

There are a lot of outputs for QPI calculation. including

arc.dat_l, arc.dat_r, arc_l.gnu, arc_r.gnu, arc.jdat_l, arc.jdat_r, arc.jsdat_l, arc.jsdat_r, arc_l_jdos.gnu, arc_l_jsdos.gnu, arc_r_jdos.gnu, arc_r_jsdos.gnu.

The gnu script with ‘only’ means we only plot the spectrum with the surface contribution, exclude the bulk contribution. (we remove file arc_l_only.gnu after v2.4.1 for the reason of misleading) jdat_l is the QPI data without consideration of spin scattering. jsdat_l is the QPI data in consideration of spin scattering.

2.3.13. Fermi arc calculation

Surface state spectrum at fixed energy level E_arc set in NAMELISTS PARAMETERS . Set SlabArc_calc=T, and set Nk1, Nk2, in NAMELISTS PARAMETERS, set k plane in KPLANE_SLAB CARD. Get the plots with “gnuplot arc_l.gnu”, “gnuplot arc_r.gnu”. _l and _r means the top and bottom surface.

2.3.14. Spin texture calculation

Spin texture calculation at a fixed energy level E_arc set in NAMELISTS PARAMETERS . Set Slabspintexture_calc=T, and set Nk1, Nk2, in NAMELISTS PARAMETERS, set k plane in KPLANE_SLAB CARD. Get the plots with “gnuplot spintext_l.gnu”, “gnuplot spintext_r.gnu”.

Note

Here we asumme that the tight-binding basis are pure spin up or pure spin down, which means that the spin up and spin down are not mixed in the basis. This could be realized if you don’t do the maximal-localized step by setting num_iter=0 in wannier90.win and select the projectors, disentanglement windown properly. If your Wannier functions are the maximal localized ones, then this feature doesn’t work. I suggest you using Wannier90 to get spin-texture which needs information from the first-principle calculations.

2.3.15. Berry phase calculation

Calculate Berry phase of a closed k path in 3D BZ. This is useful in a nodal line system. It is demonstrated that the Berry phase around a closed mirror symmetric k loop is either 0 or pi for a mirror protect nodal line system.

In WannierTools, you can specify a k path by a serials k points. Here we take the WC example, which has two nodal lines around K point.

Input

&CONTROL
BerryPhase_calc = T
/
&SYSTEM
NumOccupied = 10        ! Number of occupied Wannier orbitals
/
&PARAMETERS
Nk1 = 21    ! No. of slices for the 1st reciprocal vector
/

KPATH_BERRY
11
Direct
 0.3    0.333  -0.2
 0.3    0.333  -0.1
 0.3    0.333  -0.0
 0.3    0.333   0.1
 0.3    0.333   0.2
 0.33   0.333   0.2
 0.33   0.333   0.1
 0.33   0.333   0.0
 0.33   0.333  -0.1
 0.33   0.333  -0.2
 0.3    0.333  -0.2

Output

The value of Berry phase can be found in the WT.out.

Note

1. In principlely, the Berry phase for around a nodal line should be interger. However, the MLWF-TB model usally is not symmetric. So the value of Berry phase is close to zero or one.

2. The first and the last kpoints in the KPATH_BERRY should be the same to form a loop.

2.3.16. Berry curvature calculation for 3D bulk case

Calculate Berry curvature at a fixed k plane in 3D BZ. Set BerryCurvature_calc=T, and set Nk1, Nk2, in NAMELISTS PARAMETERS, set k plane in KPLANE_BULK CARD. Get the plot with “gnuplot Berrycurvature.gnu”.

please set NumOccpuied correctly. It represents the “occpuied” wannier bands, not the total number of electrons. In this application, the Berrycurvature is the summation over NumOccupied bands.

A typical input (take ZrTe as an example):

&CONTROL
 BerryCurvature_calc=T
/
&SYSTEM
NumOccupied = 8         ! Number of occupied Wannier orbitals
/
&PARAMETERS
Nk1 = 101    ! No. of slices for the 1st reciprocal vector
Nk2 = 101    ! No. of slices for the 2st reciprocal vector
/

KPLANE_BULK
0.00  0.00  0.00   ! Central point for 3D k slice  k3=0
1.00  0.00  0.00   ! The first vector. Integrate along this direction to get WCC
0.00  1.00  0.00   ! WCC along this direction, for Z2, usually half of the reciprocal lattice vector

2.3.17. Berry curvature calculation for slab system

Note

Not well tested.. Use it carefully.

A typical input:

&CONTROL
 BerryCurvature_slab_calc=T
/

&SYSTEM
NumOccupied = 8  ! Number of occupied Wannier orbitals of the unit cell
/

&PARAMETERS
Nk1 = 101    ! No. of slices for the 1st reciprocal vector
Nk2 = 101    ! No. of slices for the 2st reciprocal vector
/

KPLANE_SLAB
0.00  0.00         ! Central point for 3D k slice  k3=0
1.00  0.00         ! The first vector. Integrate along this direction to get WCC
0.00  1.00         ! WCC along this direction, for Z2, usually half of the reciprocal lattice vector

2.3.18. Anomalous Hall conductivity (AHC)

Calculate AHC in the energy range [OmegaMin, OmegaMax]. AHC is in unit of (Ohm*cm)^-1.

Relevant inputs are

&CONTROL
AHC_calc=T
/

&PARAMETERS
OmegaNum = 601    ! number of slices of energy
OmegaMin = -1.0   ! erergy range for AHC
OmegaMax =  1.0
Nk1 = 51   ! No. of slices for the 1st reciprocal vector
Nk2 = 51   ! No. of slices for the 2nd reciprocal vector
Nk3 = 51   ! No. of slices for the 3nd reciprocal vector
/

KCUBE_BULK
  0.00  0.00  0.00   ! Original point for 3D k plane
  1.00  0.00  0.00   ! The first vector to define 3d k space plane
  0.00  1.00  0.00   ! The second vector to define 3d k space plane
  0.00  0.00  1.00   ! The third vector to define 3d k cube

Output is sigma_ahe.txt.

2.3.19. Wannier charge center/Wilson loop calculation

Wannier charge center, which is sometimes called Wilson loop can be calculated by set WannierCenter_calc=T and set KPLANE_BULK CARD, set number of k points for two vectors is Nk1, Nk2 in NAMELISTS PARAMETERS. Notice: You should notice that the first vector in KPLANE_BULK CARD is the integration direction, this vector should be equal to one primitive reciprocal lattice vector. If you want to calculate the Z2 number, Please set the second vector to be half of the reciprocal lattice vector. You can get the Wannier charge center along the second k line. See more details In the paper written by Alexey. Soluyanov (2011). If you want to calculate the Chern number, Please set the second vector to be one primitive reciprocal lattice vector.

Note

Important: please set NumOccpuied correctly. It represents the “occpuied” wannier bands, not the total number of electrons. In this application, the Wilson loop is the trace of NumOccupied bands. If you want to study the topology between the 8th and the 9th band, then set NumOccupied=8.

Output

Outputs are wcc.dat and wcc.gnu, the format of wcc.dat is:

#         k      largestgap  sum(wcc(:,ik))      wcc(i, ik)(i=1, NumOccupied)
 0.00000000      0.60940556      0.99998388      0.00000850      0.07701431      0.07702018      0.19328973      0.19329593      0.28118760      0.28119336      0.49998615      0.50000060      0.71881052      0.71881646      0.80675987      0.80676424      0.92297767      0.92298328      0.99993530      0.99994085      0.99999935
 0.00312500      0.61256609      0.99998716      0.00030351      0.00082300      0.07688154      0.07709302      0.19117885      0.19525313      0.27952027      0.28297172      0.49188658      0.50810192      0.71703027      0.72048573      0.80480135      0.80887751      0.92290159      0.92311931      0.99908212      0.99967575
 0.00625000      0.61569946      0.99999557      0.00061525      0.00168708      0.07668131      0.07711296      0.18887026      0.19709893      0.27796336      0.28485788      0.48373617      0.51625672      0.71514220      0.72204447      0.80295562      0.81118878      0.92287871      0.92332259      0.99821913      0.99936414

 ......

The second column is the position of the largest gap of WCC. It is used for drawing a line to calculate the Z2 number (see A. Soluyanov 2011), From the fourth column to the last column, they are wcc for the occupied bands specified with “NumOccupied”. The third line is the summation of the WCC over all the “occupied” bands. It’s usefull for telling the Chern number.

Example

Take Bi2Se3 for example, we calculate the Wilson loop (WCC) at k3=0 plane, where k1, k2, k3 is in unit of reciprocal lattice vector. The you should set the particular inputs like this

&CONTROL
 WannierCenter_calc=T
/
&SYSTEM
NumOccupied = 10        ! Number of occupied Wannier orbitals
/
&PARAMETERS
Nk1 = 41    ! No. of slices for the 1st reciprocal vector
Nk2 = 41    ! No. of slices for the 2st reciprocal vector
/

KPLANE_BULK
0.00  0.00  0.00   ! Original point for 3D k slice  k3=0
1.00  0.00  0.00   ! The first vector. Integrate along this direction to get WCC
0.00  0.50  0.00   ! WCC along this direction, for Z2, usually half of the reciprocal lattice vector

For 2D materials stacked along z direction, you could think it as a 3D material with weak coupling along z direction. You can use this function to get the Z2 value at k3=0 plane to characterize the topology.

Use “gnuplot wcc.gnu” to get “wcc.eps” plot.

Here is an example.

2.3.20. Mirror Chern number calculation

At present, We can only calculate mirror Chern number for the simplest case (1. There is only one atom per atom’s type in the unit cell e.g. ZrTe. 2. kz=0 is the mirror plane we concern). For the more complex case, you can modify the source code by setting the mirror operator properly. Define your own mirror operator based on the atomic like Wannier functions in the symmetry.f90 and change the subroutine wanniercenter_mirror in wanniercenter.f90.

After properly setting of the mirror operator, you can run WannierTools with the basic parameters and the following additional parameters (Here we take ZrTe at kz=0 plane as an example)

&CONTROL
 MirrorChern_calc=T
/
&SYSTEM
NumOccupied = 8         ! Number of occupied Wannier orbitals
/
&PARAMETERS
Nk1 = 101   ! No. of slices for the 1st reciprocal vector
Nk2 = 101   ! No. of slices for the 2st reciprocal vector
/

KPLANE_BULK
0.00  0.00  0.00   ! Original point for 3D k slice  k3=0
1.00  0.00  0.00   ! The first vector. Integrate along this direction to get WCC, should be a close path
0.00  1.00  0.00   ! WCC along this direction, for Chern, usually one reciprocal lattice vector

Output

The mirror Chern number can be found in the WT.out. The WCC/Wilson loop is included in the files wcc-mirrorminus.dat and wcc-mirrorplus.dat. The gnuplot script is wcc-mirrorchernnumber.gnu. The format of wcc-mirrorplus.dat is:

#      k    sum(wcc(:,ik))      wcc(:, ik)
0.00000000      0.93401098      0.26748313      0.33122324      0.37761566      0.95768895
0.01000000      0.93458410      0.26776394      0.33149191      0.37747362      0.95785463
0.02000000      0.93515725      0.26806334      0.33205065      0.37717770      0.95786557
0.03000000      0.93572256      0.26838206      0.33288980      0.37673021      0.95772050
...

The first column is k=i/Nk2 (i=0, Nk2), we take the second vector defined in KPLANE_BULK as unit of 1. The second line is the summation of the WCC over all the “occupied/2” bands. It’s usefull for telling the Chern number. From the third column to the last column, they are wcc for the occupied/2 bands specified with “NumOccupied”.

2.3.21. Z2 number for 3D bulk materials

We can get Z2 topological index (v0, v1v2v3) from the z2 calculations of six time reversal invariant planes, i.e. (a) k1=0.0; (b) k1=0.5; (c) k2=0.0; (d) k2=0.5; (e) k3=0.0; (f) k3=0.5; Usually, you can call “Wannier charge center calculation for a plane” six times. Here we packed them up to get another function. You can set the input file like the following.

Input

The necessary tags that you should set in the wt.in

&CONTROL
Z2_3D_calc = T
/
&PARAMETERS
NumOccpuied = 18  ! No. of occupied wannier bands
Nk1 = 41    ! No. of slices of the k points for WCCs
Nk2 = 41    ! No. of slices of the k points for WCCs
/

Output

Outputs are wanniercenter3D_Z2_1.dat, wanniercenter3D_Z2_2.dat, wanniercenter3D_Z2_3.dat, wanniercenter3D_Z2_4.dat, wanniercenter3D_Z2_5.dat, wanniercenter3D_Z2_6.dat and wanniercenter3D_Z2.gnu. The z2 value can be found in the WT.out by searching “z2 number for 6 planes”. The WCC (Wilson loop) plots wanniercenter3D_Z2.eps can be got with:

gnuplot wanniercenter3D_Z2.gnu

Note

Important: please set NumOccpuied correctly. It represents the “occpuied” wannier bands, not the total number of electrons. In this application, the Wilson loop is the trace of NumOccupied bands. If you want to study the topology between the 8th and the 9th band, then set NumOccupied=8.

About the Z2 index for 3D system.

v0= (z2(ki=0)+z2(ki=0.5))mod 2

vi= z2(ki=0.5)

For the 2D system, if you set the Z axis as the stack axis, please only take the Z2 number at k3=0 plane.

2.3.22. Chern number for 3D bulk materials

Basically, you can calculate the Chern number for a closed manifold, for example, a 2D torus. For this purpose, I would suggest you using

WannierCenter_calc=T in the calculation.

We can get Chern number of six k planes, i.e. k1=0.0; k1=0.5; k2=0.0; k2=0.5; k3=0.0; k3=0.5; where k1, k2, k3 is in fractional units. Usually, you can call “Wannier charge center calculation for a plane” six times. Here we packed them up to get another function. You can set the input file like the following.

Basically, the method used here is very similar to the one used in the Z2 number calculations. We also use the Wilson loop method. However, for Z2 calculation, you only take half of the size of a time reversal invariant slice. For Chern number calculation, you need a closed momentum surface. For example, for the k1=0.0 plane, half of the time reversal invariant slice would be defined like this:

k2 is in [0, 0.5]
k3 is in [0, 1]

The full closed momentum surface would defined like this

k2 is in [0, 1]
k3 is in [0, 1]

Input

The necessary tags that you should set in the wt.in

&CONTROL
Chern_3D_calc = T
/
&PARAMETERS
NumOccpuied = 18  ! No. of occupied wannier bands
Nk1 = 41    ! No. of slices of the k points for WCCs
Nk2 = 41    ! No. of slices of the k points for WCCs
/

Output

Outputs are wanniercenter3D_Chern.dat and wanniercenter3D_Chern.gnu. The Chern number can be found in the WT.out by searching “Chern number for 6 planes”. The WCC (Wilson loop) plots wanniercenter3D_Chern.eps can be got with:

gnuplot wanniercenter3D_Chern.gnu

Note

Important: please set NumOccpuied correctly. It represents the “occpuied” wannier bands, not the total number of electrons. In this application, the Wilson loop is the trace of NumOccupied bands. If you want to study the topology between the 8th and the 9th band, then set NumOccupied=8.

For the 2D system, if you set the Z axis as the stack axis, please only take the Chern number at k3=0 plane.

2.3.23. Landau level calculations

This functionality is under testing, not released yet. Developed by QSWu and YFGuan

By applying the uniform magnetic field along one lattice vector, the Landau level spectrum can be calculated by the Peierls substitution.

Here we put one example of Graphene. The input file wt.in is like this

&TB_FILE
Hrfile = 'wannier90_hr.dat'
/

!> bulk band structure calculation flag
&CONTROL
BulkBand_calc                 = T
Hof_Butt_calc                 = T
LandauLevel_k_calc            = T
LandauLevel_wavefunction_calc = F
/

&SYSTEM
NSLAB = 200             ! the size of magnetic supercell
NumOccupied = 1         ! NumOccupied
SOC = 0                 ! soc
E_FERMI = -1.2533       ! e-fermi
/

&PARAMETERS
E_arc = 0.00        ! energy for calculate Fermi Arc
OmegaNum = 201      ! number of eigenvalues to calculate the Landau levels
Nk1 = 100           ! number k points for each line in the kpath_bulk
/

LATTICE
Angstrom
2.1377110  -1.2342080   0.0000000
0.0000000   2.4684160   0.0000000
0.0000000   0.0000000   10.000000

!> used when you want to study the projections on the orbital for each band
SELECTEDORBITALS
1  ! NumberofSelectedOrbitals without spin degeneracy
1  ! SelectedOrbitals indices without spin degeneracy

ATOM_POSITIONS
2                               ! number of atoms for projectors
Direct                          ! Direct or Cartisen coordinate
C 0.333333 0.666667 0.500000 C
C 0.666667 0.333333 0.500000 C

PROJECTORS
1 1        ! number of projectors
C  pz
C  pz

SURFACE
 0  0  1     ! magnetic field direction in units of lattice vectors
 1  0  0

KPATH_BULK            ! k point path
1              ! number of k line only for bulk band
M   0.50000  0.00000  0.00000   G   0.00000  0.00000  0.00000

WANNIER_CENTRES
Cartesian
0.712570  1.234209  5.000000
1.425141 -0.000001  5.000000

We can calculate the Hofstader butterfly plot by setting Hof_Butt_calc = T. Nslab is the size of the supercell. The corresponding magnetic field strength can be found in the output WT.out.

Since the calculation for magnetic supercell is very heavy, we have two versions of eigenvalue solvers. One is direct diagonalization, the other one is the ARPACK solver, which is based on the sparse matrix. So you have to install the ARPACK package, and specify the library in the Makefile. You don’t have to choose the solvers. It is automatically chosen according to the matrix dimensionality. If the dimensionality of the Hamiltonian matrix of the magnetic supercell is larger than 1600, WannierTools will call the sparse matrix solver.

Once the sparse matrix solver is chosen. You have to set E_arc and OmegaNum, which means you are going to calculate 2*OmegaNum+1 Landau energy levels around E_arc respect to the Fermi level.

The magnetic field is along the first vector specified in the SURFACE card.

2.3.24. Ordinary magnetoresistance calculations

After version 2.6.0, WannierTools is able to calculate magnetoresistance of a non-magnetic metal or semimetal which has Fermi surfaces. The theory part was described in Phys. Rev. B 99, 035142 (2019).

With control tag Boltz_OHE_calc=T and Symmetry_Import_calc = T, we can get the conductivity tensor for a given magnetic direction specified by Btheta and Bphi.

Here we put one example of Cu. The input file wt.in is like this

&TB_FILE
Hrfile = 'wannier90_hr.dat_nsymm48'   ! Here we symmetrized the TB model
/

&CONTROL
Boltz_OHE_calc        = T   ! calculate ordinary magnetoresistance
Symmetry_Import_calc  = T ! please set it to be true for magnetoresistance calculation
/

&SYSTEM
SOC = 0                ! There is no SOC in the hr file : SOC=0; with soc : SOC=1
E_FERMI = 7.7083       ! e-fermi in eV
Btheta= 0, Bphi= 90    ! magnetic field direction, Btheta is the angle with z axial, Bphi is the angle with respect to x axial in the x-y plane
NumOccupied = 6        ! set it anyway even don't use it. Usually, It's the valance band maximum band index.
/

&PARAMETERS
OmegaNum = 1        ! omega number
OmegaMin =  0.0     ! energy interval
OmegaMax =  0.0     ! energy interval  chemical potential \mu_i=OmegaMin+ (i-1)/(OmegaNum-1)*(OmegaMax-OmegaMin)
Nk1 =61             ! Kmesh(1) for KCUBE_BULK
Nk2 =61             ! Kmesh(2) for KCUBE_BULK
Nk3 =61             ! Kmesh(3) for KCUBE_BULK
BTauNum= 100        ! Number of B*tau we calculate
BTauMax = 40.0      ! The maximum B*tau, starting from Btau=0.
Tmin = 30           ! Temperature in Kelvin
Tmax = 330          ! Temperature in Kelvin
NumT = 11           ! number temperature we calculate. T_i=Tmin+(Tmax-Tmin)/(NumT-1)*(i-1)
/

LATTICE
Angstrom  ! Unit of length: Angstrom or Bohr
   0.0000000   1.8075000   1.8075000
   1.8075000   0.0000000   1.8075000
   1.8075000   1.8075000   0.0000000

ATOM_POSITIONS
1                               ! number of atoms for projectors
Cartisen                          ! Direct or Cartisen coordinate
Cu    0.000000      0.000000      0.000000

PROJECTORS
9            ! number of projectors
Cu s pz px py dz2 dxz dyz dx2-y2 dxy

SELECTEDBANDS
1
6   ! the 6'th band is crossing the Fermi level.

KCUBE_BULK   ! in unit of primitive reciprocal lattice vectors.
 0.00  0.00  0.00   ! Original point for 3D k plane
 1.00  0.00  0.00   ! The first vector to define 3d k space plane
 0.00  1.00  0.00   ! The second vector to define 3d k space plane
 0.00  0.00  1.00   ! The third vector to define 3d k cube

The above input file indicates that 1). The magnetic field is along z direction (Btheta=0). 2). The chemical potential is set to Zero. by OmegaNum=1; OmegaMin=Omegamax=0. 3). The temperature is set as an arrary from 30K to 330K with 30 as interval by Tmin=30; Tmax=330; NumT=11. 4). The magnetic field strenth is set as B*tau= (0, 40) with 100 points by BTauMax=40; BTauNum=100. 5). The k-mesh is set as 61*61*61.

On output, the conductivity tensor sigma/tau as a function of B*Tau for each band n, each chemical potential mu and each temperature T are presented in file sigma_band_${n}_mu_${mu}eV_T_${T}K.dat

If there are serveral band, you should write a script to add up all the conductivity tensor for all the bands, then take an inverse of it to get the resistivity tensor.

On output, there are also resistivity tensor files for each band, each chemical potential and each temperature. It’s usfull to get general information for each band.

For more illustration of example Cu, you can visit https://www.wanniertools.org/examples/ordinary-magnetoresistance-of-cu/

Note

If you obtained negative magnetoresistance at large Btau, then please increase NSlice_Btau_Max such that the negative magnetoresistance vanishes. The Lorentz force will only give us positive magnetoresistance.