Tutorial 6: Weyl point characterization in Te

In this tutorial, we calculate the chirality of a Weyl point in Te using WIEN2K density functional package. We will demonstrate possibility of characterization of Weyl points with chirality greater than 1 as discussed in Fang et al. Phys. Rev. Lett. 108, 266802 (2012) and Tsirkin et al. Phys. Rev. B 96, 045102 (2017). The entire procedure is wrapped inside the WloopPHI python code. We construct a small loop in k space (a Wilson loop) and translate the loop such that it forms a cylinder which encloses the Weyl point of interest. WloopPHI creates a series of Wilson loops along the translation trajectory and calls BerryPI to calculate the Berry phase along each Wilson loop. More details about the method can be found in Comp. Phys. Commun. 270, 108147 (2022). We are thankful to Stepan Tsirkin for helpful discussions and suggestions during preparation of this tutorial.

Step-by-step instructions:

1. Construct a structure file using w2web (or a makestruct utility in the command line):

   1 atom Te(0.2633, 0, 1/3), space group 152, lattice parameters (in Bohrs): a = b = 8.422887, c = 11.204186, and angles α = β = 90°, γ = 120°.

   

   Once the structure is saved, two additional equivalent Te atoms will be added automatically at (0, 0.2633, 2/3) and (0.7367, 0.7367, 0).

2. Initialize the calculation with PBE, 3% reduction of muffin-tin radii R_MT, the product min(R_MT)K_max = 7, and 500 k points in the whole BZ.

   init_lapw -b -vxc 13 -red 3 -rkmax 7 -numk 500 -nohdlo

   The -nohdlo switch is needed to avoid generation of high-derivative local orbitals by LSTART. Wien2wannier is incompatible with high-derivative local orbitals.

3. Initialize the spin-orbit coupling (SOC) using all the default parameters.

   init_so_lapw

4. Perform a self-consistent field (SCF) calculation with the energy convergence of 10⁻⁴ Ry, and charge convergence of 10⁻³ e.

   run_lapw -ec 0.0001 -cc 0.001 -so

5. Generate a A-Γ−A k path with 100 intermediate points using XCrysDen, and save the k point list to a file case.klist_band.

   

6. Recalculate the eigenvalues on the k path selected:

   x lapw1 -band
   x lapwso
   

7. In w2web under Tasks > Bandstructure edit case.insp. Insert the correct Fermi energy and select the energy range for plotting between -5 and -2 eV, set major ticks to 1 eV. Execute x spaghetti with -so option to include SOC and then plot bandstructure.

   

   The band crossing of interest is marked as WP on the figure above. We can compare it to Fig. 5a in Tsirkin et al. Phys. Rev. B 96, 045102 (2017), where each band crossing is marked with its own chirality. The one we selected corresponds to chirality of 3.

8. Determine index of the band located immediately below the crossing of interest in the figure above. For this we need to check eigenvalues at Γ point in the case.energyso file and take into account the energy range of bands listed in case.scf2 file. Note that energies on the band structure plot are measured relative to the Fermi energy.

9. We need to construct a Wilson loop as shown in the figure below.

   

   This python utility WloopIN.py is very helpful. We will need to edit this section to set up the Weyl point tentative location, loop radius, loop translation direction, starting and ending coordinates along the translation direction, number of points on the loop:

   # Enter the precise location of Weyl point (x, y, z).
   cx, cy, cz = 0, 0, 0
   # Enter the radius of a loop
   r = 0.1
   # Trajectory direction.
   direction = str('z')
   # Starting point and End point.
   S, E = -0.15, 0.15
   # How many vertex points for Wilson loop
   n = 16
   

   The script will prepare a file WloopCoordinate.dat with loop coordinates.

10. In the case folder create the input file Wloop.in. It specifies the loop coordinates (prepared in the previous step), number of steps along the translation direction, and the band range.

    31
    1:39
    &WloopCoordinate
    0.10000 0.00000 -0.15000 ; 0.10000 0.00000 0.15000
    0.09239 0.03827 -0.15000 ; 0.09239 0.03827 0.15000
    0.07071 0.07071 -0.15000 ; 0.07071 0.07071 0.15000
    0.03827 0.09239 -0.15000 ; 0.03827 0.09239 0.15000
    0.00000 0.10000 -0.15000 ; 0.00000 0.10000 0.15000
    -0.03827 0.09239 -0.15000 ; -0.03827 0.09239 0.15000
    -0.07071 0.07071 -0.15000 ; -0.07071 0.07071 0.15000
    -0.09239 0.03827 -0.15000 ; -0.09239 0.03827 0.15000
    -0.10000 0.00000 -0.15000 ; -0.10000 0.00000 0.15000
    -0.09239 -0.03827 -0.15000 ; -0.09239 -0.03827 0.15000
    -0.07071 -0.07071 -0.15000 ; -0.07071 -0.07071 0.15000
    -0.03827 -0.09239 -0.15000 ; -0.03827 -0.09239 0.15000
    -0.00000 -0.10000 -0.15000 ; -0.00000 -0.10000 0.15000
    0.03827 -0.09239 -0.15000 ; 0.03827 -0.09239 0.15000
    0.07071 -0.07071 -0.15000 ; 0.07071 -0.07071 0.15000
    0.09239 -0.03827 -0.15000 ; 0.09239 -0.03827 0.15000
    END
    

    BerryPI will automatically connect the first and last points to form a loop in k space. It is tricky to select a correct size of the loop. The loop should be small enough to capture only one Weyl point. However, tiny loops require more steps along the translation direction to achieve a smooth phase evolution.

11. Execute WloopPHI from the WIEN2k case directory:

    python $WIENROOT/SRC_BerryPI/BerryPI/WloopPHI.py Wloop.in

    This is the most computationally intensive part. The WloopPHI script automatically executes the BerryPI to calculate the Berry phase φ on each Wilson loop along the translation trajectory. The output data are written in the PHI.dat file.

    # Loop (z)       BerryPhase(BP)   BP(-/+pi wrap)   BP(unwrap)/2pi
    -0.15000          -7.23112          -0.94794          -0.15087
    -0.14000          5.20487          -1.07831          -0.17162
    -0.13000          11.32915          -1.23722          -0.19691
    -0.12000          4.85522          -1.42797          -0.22727
    -0.11000          -1.65461          -1.65461          -0.26334
    -0.10000          -8.20733          -1.92415          -0.30624
    -0.09000          -2.24443          -2.24443          -0.35721
    -0.08000          3.65709          -2.62609          -0.41796
    -0.07000          -9.36478          -3.08159          -0.49045
    ...
    

12. Plot the Berry phase vs k_z from PHI.dat file (first and last columns). The Weyl point chirality is the phase (divided by 2π) accumulated while the loop traverses through the Weyl point. It is the same as a flux of a Berry curvature through the wall of a cylinder that encloses the Weyl point. The flux through a closed surface that contains a Weyl point inside can only be an integer number multiple of 2π.

    

    The initial phase is not exactly 0 since the loop starts and finishes close the the Weyl point (k_z=-0.15 and 0.15) when the flux of the Berry curvature through the loop is still finite. We could have started further away, but it would require more steps along the translation direction to achieve a smooth evolution of the Berry phase vs k_z.

    

    If we count the flux through all surfaces of the cylinder, we can conclude that the Weyl point chirality is = 3 which is consistent with Tsirkin et al. Phys. Rev. B 96, 045102 (2017). It is "unusual" for Weyl points to have the chirality different from 1 (or -1). Our band crossing of interest has a linear dispersion in the Gamma-A direction and a unique (cubic) band dispersion in the transverse direction (Gamma-M).

    

    The sign of the phase is a matter of convention. In this case our phase winds down. If we encounter a Weyl point of the opposite chirality, the phase will wind up provided we are consistent with the order/direction in which we select k points on the loop and the translation direction. (See discussion of the right-hand rule in Comp. Phys. Commun. 270, 108147 (2022).)