3. Band structure

In this tutorial, we will use Wananier interpolation to calculate the band structure of silicon valence and conduction manifold. A bit different from previous tutorials, we will

construct a InterpModel by reading the win, mmn, eig, and chk.fmt files
run interpolate on the InterpModel to compute band structure
read Wannier90 interpolated band.dat and compare with Wannier.jl interpolated bands

Tip

This is a HTML version of the tutorial, you can download corresponding

Jupyter notebook: tutorial.ipynb
Julia script: tutorial.jl

Preparation

Load the package

using Wannier
using WannierPlots

The WannierPlots is a companion package of Wannier.jl, for plotting band structures, real space WFs, etc. This separation reduces the compiling latency of Wannier.jl, and avoids introduing too many package dependencies that are irrelevant to Wannierization.

Tip

Use the run.sh script which automate the scf, nscf, pw2wannier90 steps.

Model generation

from `chk.fmt` file

We will use the read_w90_interp function to read the win, mmn, eig, and chk.fmt files, and construct an InterpModel that is used for interpolation purpose.

model = read_w90_interp("si2")


recip_lattice: Å⁻¹
  a1: -1.15701  1.15701  1.15701
  a2:  1.15701 -1.15701  1.15701
  a3:  1.15701  1.15701 -1.15701
kgrid  : 4 4 4
n_kpts : 64

lattice: Å
  a1:  0.00000  2.71527  2.71527
  a2:  2.71527  0.00000  2.71527
  a3:  2.71527  2.71527  0.00000
grid    = 4 4 4
n_rvecs = 93
using MDRS interpolation

KPath{3} (6 points, 1 paths, 8 points in paths):
 points: :U => [0.625, 0.25, 0.625]
         :G => [0.0, 0.0, 0.0]
         :W => [0.5, 0.25, 0.75]
         :K => [0.375, 0.375, 0.75]
         :L => [0.5, 0.5, 0.5]
         :X => [0.5, 0.0, 0.5]
  paths: [:G, :X, :U, :K, :G, :L, :W, :X]
  basis: [-1.157011, 1.157011, 1.157011]
         [1.157011, -1.157011, 1.157011]
         [1.157011, 1.157011, -1.157011]

Tip

To avoid a bloated Model that contains everything, and to "Do One Thing And Do It Well", we separate on purpose the Model that is solely for Wannierization, and the InterpModel, that is only used for Wannier interpolation of operators. This is convenient for developers to focus on the the Wannierization or interpolation algorithm without being distracted by the other part.

from `amn` file

You can also use the read_w90 function to read the amn file and generate a Model,

m = read_w90("si2")

lattice: Å
  a1:  0.00000  2.71527  2.71527
  a2:  2.71527  0.00000  2.71527
  a3:  2.71527  2.71527  0.00000

atoms: fractional
  Si:  0.00000  0.00000  0.00000
  Si:  0.25000  0.25000  0.25000

n_bands: 16
n_wann : 8
kgrid  : 4 4 4
n_kpts : 64
n_bvecs: 8

b-vectors:
         [bx, by, bz] / Å⁻¹                weight
  1       0.28925   -0.28925    0.28925    1.49401
  2      -0.28925   -0.28925   -0.28925    1.49401
  3      -0.28925    0.28925   -0.28925    1.49401
  4      -0.28925   -0.28925    0.28925    1.49401
  5       0.28925    0.28925   -0.28925    1.49401
  6      -0.28925    0.28925    0.28925    1.49401
  7       0.28925   -0.28925   -0.28925    1.49401
  8       0.28925    0.28925    0.28925    1.49401

Then use the InterpModel constructor to construct an InterpModel from an existing Model,

m = Wannier.InterpModel(m)


recip_lattice: Å⁻¹
  a1: -1.15701  1.15701  1.15701
  a2:  1.15701 -1.15701  1.15701
  a3:  1.15701  1.15701 -1.15701
kgrid  : 4 4 4
n_kpts : 64

lattice: Å
  a1:  0.00000  2.71527  2.71527
  a2:  2.71527  0.00000  2.71527
  a3:  2.71527  2.71527  0.00000
grid    = 4 4 4
n_rvecs = 93
using MDRS interpolation

KPath{3} (6 points, 2 paths, 8 points in paths):
 points: :U => [0.625, 0.25, 0.625]
         :W => [0.5, 0.25, 0.75]
         :K => [0.375, 0.375, 0.75]
         :Γ => [0.0, 0.0, 0.0]
         :L => [0.5, 0.5, 0.5]
         :X => [0.5, 0.0, 0.5]
  paths: [:Γ, :X, :U]
         [:K, :Γ, :L, :W, :X]
  basis: [-1.157011, 1.157011, 1.157011]
         [1.157011, -1.157011, 1.157011]
         [1.157011, 1.157011, -1.157011]

However, there are some differences:

the amn gauge is used, instead of that from chk
the kpoint path for band structure is auto generated from the lattice, instead of using that in win file (if found)

So, it is recommended to use the read_w90_interp function, or you run a max_localize or disentangle on the Model to smooth the gauge, then construct an InterpModel.

Tips on kpath

The read_w90_interp function will parse the kpoint_path block in the win file:

if the kpoint_path block is found, it will use that path
if the kpoint_path block is not found, the InterpModel constructor will auto generate a kpath from the lattice, by using the get_kpath function

Tip

The get_kpath works for arbitrary lattice, e.g., either conventional or primitive, it will generate the correct kpath.

During the band interpolation, the interpolate function will call the interpolate_w90 function to generate an exactly identical kpath to that of Wannier90.

Tips on interpolation algorithm

There are two interpolation algorithms

Wigner-Seitz (WS) interpolation
Minimal-distance replica selection (MDRS) method

The read_w90_interp function will parse the use_ws_distance parameter in the win file:

if use_ws_distance is true, it will use the MDRS
if use_ws_distance is false, it will use the WS
if not found, it will use the MDRS

You can also control the interpolation algorithm by the constructor InterpModel(model::Model; mdrs::Bool=true).

Moreover, there are two versions of MDRS, i.e., mdrsv1 and mdrsv2 in the jargon of Wannier.jl. The v2 is faster than v1 so it is used by default. For details, see the API for interpolation, e.g., fourier, invfourier, _fourier_mdrs_v1, _fourier_mdrs_v2, _invfourier_mdrs_v1, and _invfourier_mdrs_v2, etc.

If you are tech-savvy, you can even control the generation of R vectors for interpolation, this is useful for developing new interpolation algorithms, see get_Rvectors_mdrs, get_Rvectors_ws, etc.

Band interpolation

Computing band structure is very easy, by calling the interpolate function, it will return a KPathInterpolant, and the band eigen energies,

kpi, E = interpolate(model)

(StaticArraysCore.SVector{3, Float64}[[0.0, 0.0, 0.0], [0.005, 0.0, 0.005], [0.01, 0.0, 0.01], [0.015, 0.0, 0.015], [0.02, 0.0, 0.02], [0.024999999999999998, 0.0, 0.024999999999999998], [0.03, 0.0, 0.03], [0.034999999999999996, 0.0, 0.034999999999999996], [0.04, 0.0, 0.04], [0.04499999999999999, 0.0, 0.04499999999999999]  …  [0.5, 0.045000000000000005, 0.545], [0.5, 0.04, 0.5399999999999999], [0.5, 0.034999999999999996, 0.535], [0.5, 0.029999999999999985, 0.53], [0.5, 0.02499999999999998, 0.525], [0.5, 0.020000000000000025, 0.52], [0.5, 0.015000000000000017, 0.515], [0.5, 0.010000000000000012, 0.51], [0.5, 0.005000000000000006, 0.505], [0.5, 0.0, 0.5]], [-5.798515544176985 -5.797983031391298 … -1.6351163088568883 -1.6353082316090033; 6.202325333824968 6.198803664460466 … -1.635116304066449 -1.6353082316089598; … ; 8.70347534455524 8.706030783974015 … 16.246232171720404 16.247897361941128; 9.783552118546163 9.786953281936183 … 16.24623220558102 16.247897361941433])

The kpi stores the specific kpoint coordinates along the kpath, while the InterpModel.kpath only stores the high-symmetry kpoints their labels, that's why we need to return the kpi object

kpi

511-element Brillouin.KPaths.KPathInterpolant{3}:
 [0.0, 0.0, 0.0]
 [0.005, 0.0, 0.005]
 [0.01, 0.0, 0.01]
 [0.015, 0.0, 0.015]
 [0.02, 0.0, 0.02]
 [0.024999999999999998, 0.0, 0.024999999999999998]
 [0.03, 0.0, 0.03]
 [0.034999999999999996, 0.0, 0.034999999999999996]
 [0.04, 0.0, 0.04]
 [0.04499999999999999, 0.0, 0.04499999999999999]
 ⋮
 [0.5, 0.04, 0.5399999999999999]
 [0.5, 0.034999999999999996, 0.535]
 [0.5, 0.029999999999999985, 0.53]
 [0.5, 0.02499999999999998, 0.525]
 [0.5, 0.020000000000000025, 0.52]
 [0.5, 0.015000000000000017, 0.515]
 [0.5, 0.010000000000000012, 0.51]
 [0.5, 0.005000000000000006, 0.505]
 [0.5, 0.0, 0.5]

Plotting band structure

Saving interpolated band

You can save the result to the same format as Wannier90 band.dat, by

write_w90_band("wjl", kpi, E)

[ Info: Written to wjl_band.kpt
[ Info: Written to wjl_band.dat
[ Info: Written to wjl_band.labelinfo.dat

where wjl is the seedname of the output, i.e., written files are

wjl_band.kpt
wjl_band.dat
wjl_band.labelinfo.dat

Visualization in the Julia world

Instead of saving, you can also plot the band structure by calling the WannierPlots.plot_band function,

Note

This requires you having executed using PlotlyJS

P = plot_band(kpi, E)

Or, you can use the plotting functions provided by Brillouin.jl, but requires a little bit transposation to size (n_kpts, n_bands)

P = plot(kpi, E')

Comparing band structures

Now we load the Wannier90 interpolated band, to compare between the two codes,

kpi_w90, E_w90 = read_w90_band("si2", model.recip_lattice)

(StaticArraysCore.SVector{3, Float64}[[0.0, 0.0, 0.0], [0.005, 0.0, 0.005], [0.01, 0.0, 0.01], [0.015, 0.0, 0.015], [0.02, 0.0, 0.02], [0.025, 0.0, 0.025], [0.03, 0.0, 0.03], [0.035, 0.0, 0.035], [0.04, 0.0, 0.04], [0.045, 0.0, 0.045]  …  [0.5, 0.045, 0.545], [0.5, 0.04, 0.54], [0.5, 0.035, 0.535], [0.5, 0.03, 0.53], [0.5, 0.025, 0.525], [0.5, 0.02, 0.52], [0.5, 0.015, 0.515], [0.5, 0.01, 0.51], [0.5, 0.005, 0.505], [0.5, 0.0, 0.5]], [-5.7985155 -5.797983 … -1.6351163 -1.6353082; 6.2023253 6.1988037 … -1.6351163 -1.6353082; … ; 8.7034753 8.7060308 … 16.246232 16.247897; 9.7835521 9.7869533 … 16.246232 16.247897])

Tip

Here I pass a recip_lattice to the read_w90_band function, so that it will return a tuple of (KPathInterpolant, Matrix). You can also call the read_w90_band function without recip_lattice, however, this "raw" version will return rather verbose outputs, not very handy for usage. See the API read_w90_band for details.

Now let's have a look at Wannier90 interpolated band

P = plot_band(kpi_w90, E_w90)

Finally, do the comparison

P = plot_band_diff(kpi, E_w90, E)

As expected, the two band structures exactly overlaps. 🥳

Well done! You have finished the first Wannier interpolation.

This page was generated using Literate.jl.