9. Parallel transport using custom b-vectors

In this tutorial, we will Wannierize a single isolated band of CuBr2, by constructing the parallel transport gauge (PTG).

However, the CuBr2 system has a special set of b-vectors: some nearest neighbors are not included. This breaks the requirement of the PTG, since it needs the overlap matrices of nearest neighbors along 6 directions. Therefore, we need to generate a custom set of b-vectors, write them to a nnkp file, rerun pw2wannier90.x, and use the new mmn file to construct the PTG.

Warning

The generated b-vectors do not satisfy the completeness condition, therefore it should not be used to compute the spread and other physical operators.

Outline

  1. construct a Model, by reading the win, amn, mmn, and eig files
  2. generate a new set of b-vectors containing only the 6 directions
  3. write the new b-vectors to a nnkp file
  4. launch pw2wannier90.x to compute a new mmn file
  5. read the new mmn file, and construct the PTG
  6. use the PTG gauge matrices in the Model
  7. maxiaml localize to smoothen the gauge
  8. interpolate band structure
Tip

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

Preparation

Load the package

using Wannier
using WannierPlots

Model generation

We will use the read_w90 function to read the win, amn, mmn, and eig files, and construct a Model

Tip

Before you do the actual calculations, if you don't know the kpath of your structure, you can use the following functions to auto generate a kpath

win = read_win("CuBr2.win")
kp = get_kpath(win.unit_cell, win.atoms_frac, win.atom_labels)
kpi = Wannier.interpolate_w90(kp, 100)
Wannier.write_w90_kpt_label("CuBr2", kpi)

The outputs are two files: CuBr2_band.kpt and CuBr2_band.labelinfo.dat, which are exactly in the same format as the Wannier90 output bands, you can inspect these two files and use them to populate the kpoint_path block in the win file, or the K_POINTS block in the pw.x bands calculation.

model = read_w90("CuBr2")
lattice: Å
  a1:  7.14000  0.00000  0.00000
  a2:  1.85704  1.73000 -3.07238
  a3: -1.85704  1.73000  3.07238

atoms: fractional
  Cu:  0.00000  0.50000  0.50000
  Br:  0.24000  0.98500  0.01500
  Br:  0.76000  0.01500  0.98500

n_bands: 1
n_wann : 1
kgrid  : 4 4 4
n_kpts : 64
n_bvecs: 10

b-vectors:
         [bx, by, bz] / Å⁻¹                weight
  1       0.22000    0.00000    0.13297    7.49333
  2      -0.22000    0.00000   -0.13297    7.49333
  3       0.22000    0.00000   -0.37829    2.53647
  4      -0.22000    0.00000    0.37829    2.53647
  5       0.44000    0.00000   -0.24532   -0.53128
  6      -0.44000    0.00000    0.24532   -0.53128
  7      -0.22000    0.45399    0.12266    1.21298
  8       0.22000   -0.45399   -0.12266    1.21298
  9       0.22000    0.45399   -0.12266    1.21298
 10      -0.22000   -0.45399    0.12266    1.21298
Note

This is a simple system, you can even get good WFs with a simple max_localize(model). However, we will use this tutorial to demonstrate the method to overcome the b-vector problem in PTG.

The initial spread is

omega(model)
  WF     center [rx, ry, rz]/Š             spread/Ų
   1     1.30125    -1.18938     0.77896    95.06736
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
   ΩI  =     5.16812
   Ω̃   =    89.89923
   ΩOD =     0.00000
   ΩD  =    89.89923
   Ω   =    95.06736

Custom b-vectors

As shown here, the b-vectors are a bit strange that some nearest neighbors are not included,

using PlotlyJS
WannierPlots.plot(model.bvectors)

thus we need to use a different set of b-vectors, which simply includes all the 6 nearest neighbors.

bvectors_nn = Wannier.get_bvectors_nearest(model.kpoints, model.recip_lattice)
WannierPlots.plot(bvectors_nn)

now write to a new nnkp file

n_wann = model.n_wann
exclude_bands = collect(1:16)
Wannier.write_nnkp("CuBr2_nn.nnkp", bvectors_nn, n_wann, exclude_bands)
[ Info: Writing nnkp file: CuBr2_nn.nnkp

now we can run pw2wannier90.x to generate a new mmn file. A template input file for pw2wannier90.x is provided here

&INPUTPP
    outdir = './out/'
    prefix = 'cubr2'
    seedname = 'CuBr2_nn'
    write_amn = .false.
    write_eig = .false.
    scdm_proj = .true.  ! this is useless, but we need it to avoid p2w error
/

Parallel transport

Now, we load the new mmn file, and construct a new Model

M_nn = WannierIO.read_mmn("CuBr2_nn.mmn")[1];
model_nn = Wannier.Model(
    model.lattice,
    model.atom_positions,
    model.atom_labels,
    model.kgrid,
    model.kpoints,
    bvectors_nn,
    model.frozen_bands,
    M_nn,
    model.U,
    model.E,
)
lattice: Å
  a1:  7.14000  0.00000  0.00000
  a2:  1.85704  1.73000 -3.07238
  a3: -1.85704  1.73000  3.07238

atoms: fractional
  Cu:  0.00000  0.50000  0.50000
  Br:  0.24000  0.98500  0.01500
  Br:  0.76000  0.01500  0.98500

n_bands: 1
n_wann : 1
kgrid  : 4 4 4
n_kpts : 64
n_bvecs: 6

b-vectors:
         [bx, by, bz] / Å⁻¹                weight
  1       0.22000    0.00000    0.13297    1.00000
  2      -0.22000    0.00000   -0.13297    1.00000
  3       0.00000    0.45399   -0.25563    1.00000
  4       0.00000   -0.45399    0.25563    1.00000
  5       0.00000    0.45399    0.25563    1.00000
  6       0.00000   -0.45399   -0.25563    1.00000

finally, we run parallel transport

U, _ = parallel_transport(model_nn)
([1.0 - 7.705582489488652e-17im;;; 0.35424192344475075 + 0.935153816050688im;;; 0.8320062948033026 - 0.5547661898563032im;;; … ;;; 0.6944757027862014 - 0.7195161556487887im;;; -0.06397879638386897 - 0.9979512581350206im;;; -0.9846154177322565 - 0.17473545479934413im], Wannier.Obstruction{ComplexF64}([1.0000000000000002 - 5.551115123125783e-17im;;; 1.0000000000000004 + 0.0im;;; 0.9999999999999998 - 1.1102230246251565e-16im;;; 1.0000000000000002 - 5.551115123125783e-17im], [0.7071068004161116 + 0.707106761956983im;;; 0.7071066324925186 + 0.7071069298805451im;;; 0.7071067413965386 + 0.7071068209765542im;;; 0.7071065342094165 + 0.7071070281635922im], [0.980785281729587 + 0.19509031534808527im;;; 0.98078529056235 + 0.19509027094277775im;;; 0.9807852760386562 + 0.19509034395832411im;;; 0.9807852791050835 + 0.19509032854235353im], [1.0 + 0.0im;;; 1.0 + 0.0im;;; 1.0 + 0.0im;;; 1.0 + 0.0im;;;; 1.0 - 1.3877787807814457e-17im;;; 1.0 + 0.0im;;; 1.0 - 2.7755575615628926e-17im;;; 1.0 - 1.3877787807814457e-17im;;;; 1.0 - 2.7755575615628907e-17im;;; 1.0 + 0.0im;;; 1.0 - 5.551115123125783e-17im;;; 1.0 - 2.7755575615628907e-17im;;;; 1.0 - 4.163336342344336e-17im;;; 1.0 + 0.0im;;; 1.0 - 8.326672684688675e-17im;;; 1.0 - 4.163336342344336e-17im], [0.9999999999999999 - 8.854311322427737e-17im;;; 0.9999999999999999 - 1.6704774771121138e-16im;;; 1.0 + 2.2479189238238287e-17im;;; 0.9999999999999999 - 1.6704774771121138e-16im;;;; 0.9807852817295866 + 0.1950903153480853im;;; 0.9807852701470825 + 0.1950903735772582im;;; 0.9807852776587205 + 0.19509033581370977im;;; 0.9807852633680203 + 0.19509040765789296im;;;; 0.9238795377147694 + 0.3826834198027715im;;; 0.9238794922749718 + 0.38268352950413315im;;; 0.923879521744188 + 0.38268345835916495im;;; 0.9238794656797542 + 0.3826835937106419im;;;; 0.8314696236339739 + 0.5555702160609205im;;; 0.8314695246814325 + 0.5555703641538422im;;; 0.8314695888554381 + 0.5555702681106762im;;; 0.8314694667660272 + 0.5555704508303322im], [1.0000000000000002 + 1.1487288329390876e-17im;;; 1.0 + 5.39108863965794e-18im;;; 1.0000000000000002 + 1.1487288329390876e-17im;;; 1.0 + 5.39108863965794e-18im;;;; 0.9987954562885715 + 0.04906767262979589im;;; 0.9987954568439591 + 0.04906766132461526im;;; 0.9987954559307362 + 0.04906767991369822im;;; 0.9987954561235474 + 0.04906767598893461im;;;; 0.9951847270053904 + 0.09801713694659044im;;; 0.9951847292242656 + 0.0980171144179675im;;; 0.9951847255757738 + 0.09801715146173123im;;; 0.9951847263460898 + 0.09801714364058078im;;;; 0.989176510712961 + 0.1467304694115425im;;; 0.9891765156954052 + 0.14673043582262568im;;; 0.9891765075027827 + 0.14673049105280553im;;; 0.9891765092325135 + 0.14673047939190764im]))

Since the b-vectors are not complete, computing spread on top of model_nn is meaningless, we assign back the gauge matrices to model

model.U .= U;

and the new spread is

omega(model)
  WF     center [rx, ry, rz]/Š             spread/Ų
   1    -0.00000    -1.73000     0.00000     6.63661
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
   ΩI  =     5.16812
   Ω̃   =     1.46849
   ΩOD =     0.00000
   ΩD  =     1.46849
   Ω   =     6.63661

then further maximal localize it

U = max_localize(model)
model.U .= U;
omega(model)
  WF     center [rx, ry, rz]/Š             spread/Ų
   1    -0.00000    -1.73000     0.00000     5.16812
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
   ΩI  =     5.16812
   Ω̃   =     0.00000
   ΩOD =     0.00000
   ΩD  =     0.00000
   Ω   =     5.16812

Band interpolation

Finally, let's have a look at the band interpolation.

First load the QE band structure

kpoints_qe, E_qe = WannierIO.read_qe_band("qe_bands.dat");

the Fermi energy from scf calculation

ef = 4.6459
4.6459

Force using kpoint_path in win file

win = read_win("CuBr2.win")
kpath = Wannier.get_kpath(win.unit_cell, win.kpoint_path)

interp_model = Wannier.InterpModel(model; kpath=kpath)

interpolate band structure the QE bands use 50 points per segment, so we use 50 here as well

kpi = Wannier.interpolate_w90(kpath, 50)
E = Wannier.interpolate(interp_model, kpi)

plot band difference

fig = plot_band_diff(kpi, E_qe, E; fermi_energy=ef)
fig.plot.layout.width = 500
fig.plot.layout.height = 500
fig.plot.layout.autosize = false
fig
Main.HTMLPlot(fig, 500)  # hide

That's all!

This page was generated using Literate.jl.