1. Maximal localization of isolated manifold

In the first tutorial, we run the maximal localization algorithm on the silicon valence bands. We need to

generate the amn, mmn, and eig files by using Quantum ESPRESSO (QE)
construct a Model for Wannier.jl, by reading the win, amn, mmn, and eig files
run Wannier.jl max_localize on the Model to minimize the spread
write the maximal localized gauge to a new amn file

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 Printf  # for pretty print

Before we running the calculation, one thing worth mentioning is the initial projection block in the win file for valence bands of silicon. In such cases, we cannot use the s and p projections, but using bond-centered s orbitals as guides for the bonding WFs.

The bond centers can be calculated by using the find_nearests function in Wannier.jl

win = read_win("si2.win")

(wannier_plot_supercell = 3, num_bands = 4, num_iter = 4000, wannier_plot = true, num_cg_steps = 200, conv_window = 3, atom_labels = ["Si", "Si"], num_wann = 4, conv_tol = 2.0e-7, wvfn_formatted = true, atoms_frac = [0.0 0.25; 0.0 0.25; 0.0 0.25], kpoints = [0.0 0.0 … 0.75 0.75; 0.0 0.0 … 0.75 0.75; 0.0 0.25 … 0.5 0.75], mp_grid = [4, 4, 4], unit_cell_cart = [0.0 2.715265 2.715265; 2.715265 0.0 2.715265; 2.715265 2.715265 0.0], spn_formatted = true, kpoint_path = Vector{Pair{Symbol, StaticArraysCore.SVector{3, Float64}}}[[:G => [0.0, 0.0, 0.0], :X => [0.5, 0.0, 0.5]], [:X => [0.5, 0.0, 0.5], :U => [0.625, 0.25, 0.625]], [:U => [0.625, 0.25, 0.625], :K => [0.375, 0.375, 0.75]], [:K => [0.375, 0.375, 0.75], :G => [0.0, 0.0, 0.0]], [:G => [0.0, 0.0, 0.0], :L => [0.5, 0.5, 0.5]], [:L => [0.5, 0.5, 0.5], :W => [0.5, 0.25, 0.75]], [:W => [0.5, 0.25, 0.75], :X => [0.5, 0.0, 0.5]]], projections = ["c=0.67882, -0.67882, -0.67882:s", "c=-0.67882, -0.67882, 0.67882:s", "c=-0.67882, 0.67882, -0.67882:s", "c=0.67882, 0.67882, 0.67882:s"], bands_plot = true, fermi_energy = 6.5589)

Tip

Refer to the API documentation for more details about the function usage.

We will find the 4 nearest atoms to the 1st atom, i.e. 4 of the periodic images of the 2nd atom, then calculate the middle point between each of these 4 atoms and the 1st atom.

Note

In Julia, array index starts from 1, and array is column-major

fractional coordiantes of the 1st Si atom

atom1 = win.atoms_frac[:, 1]

3-element Vector{Float64}:
 0.0
 0.0
 0.0

cartesian coordiantes of atom1

atom1_cart = win.unit_cell * atom1

find 4 nearest neighbors of atom1, I use 5 because the 1st returned neighbor is the atom1 itself

distances, indexes, translations = Wannier.find_nearests(
    atom1, 5, win.unit_cell, win.atoms_frac
)

print the nearest atom and bond center, in cartesian coordinates

for i in 2:5
    idx = indexes[i]
    nn_frac = win.atoms_frac[:, idx] + translations[:, i]
    nn_cart = win.unit_cell * nn_frac
    @printf(
        "nearest atom: \n  index = %d, distance = %.5f, (x, y, z) = (%.5f, %.5f, %.5f)\n",
        idx,
        distances[i],
        nn_cart...
    )
    c = (atom1_cart + nn_cart) / 2
    @printf("  bond center = (%.5f, %.5f, %.5f)\n", c...)
end

Now we can use these outputs for writing projection block in win file, and run the QE calculations for amn, mmn, and eig files.

Tip

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

Model generation

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

model = 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: 4
n_wann : 4
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

Maximal localization

Maximal localization can be easily done by calling the max_localize function, which returns the gauge matrices U,

U = max_localize(model);

[ Info: Initial spread
  WF     center [rx, ry, rz]/Å              spread/Å²
   1     0.67882    -0.67882    -0.67882     1.62976
   2    -0.67882    -0.67882     0.67882     1.62976
   3    -0.67882     0.67882    -0.67882     1.62976
   4     0.67882     0.67882     0.67882     1.62976
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
   ΩI  =     5.94279
   Ω̃   =     0.57624
   ΩOD =     0.57624
   ΩD  =     0.00000
   Ω   =     6.51903


Iter     Function value   Gradient norm
     0     6.519031e+00     7.283433e-01
 * time: 0.0012469291687011719
     1     6.517641e+00     9.528811e-04
 * time: 0.8107109069824219
     2     6.517622e+00     9.942617e-05
 * time: 0.8176000118255615
     3     6.517621e+00     5.505746e-06
 * time: 0.827380895614624
[ Info: Final spread
  WF     center [rx, ry, rz]/Å              spread/Å²
   1     0.67882    -0.67882    -0.67882     1.62941
   2    -0.67882    -0.67882     0.67882     1.62941
   3    -0.67882     0.67882    -0.67882     1.62941
   4     0.67882     0.67882     0.67882     1.62941
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
   ΩI  =     5.94279
   Ω̃   =     0.57483
   ΩOD =     0.57483
   ΩD  =     0.00000
   Ω   =     6.51762

The initial spread is

omega(model)

  WF     center [rx, ry, rz]/Å              spread/Å²
   1     0.67882    -0.67882    -0.67882     1.62976
   2    -0.67882    -0.67882     0.67882     1.62976
   3    -0.67882     0.67882    -0.67882     1.62976
   4     0.67882     0.67882     0.67882     1.62976
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
   ΩI  =     5.94279
   Ω̃   =     0.57624
   ΩOD =     0.57624
   ΩD  =     0.00000
   Ω   =     6.51903

The final spread is

omega(model, U)

  WF     center [rx, ry, rz]/Å              spread/Å²
   1     0.67882    -0.67882    -0.67882     1.62941
   2    -0.67882    -0.67882     0.67882     1.62941
   3    -0.67882     0.67882    -0.67882     1.62941
   4     0.67882     0.67882     0.67882     1.62941
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
   ΩI  =     5.94279
   Ω̃   =     0.57483
   ΩOD =     0.57483
   ΩD  =     0.00000
   Ω   =     6.51762

Since we have very good initial guess, the spread only decreases a little bit.

Note

The convergence thresholds is determined by the keyword arguments of max_localize, e.g., f_tol for the tolerance on spread, and g_tol for the tolerance on the norm of spread gradient, etc. You can use stricter thresholds to further minimize a bit the spread.

Save the new gauge

We can further do band interpolation in Wannier.jl with the new U matrices, however, for this very first tutorial, we will just save the new gauge to an amn file, which can be used as the new initial guess for Wannier90, or reuse it in Wannier.jl.

write_amn("si2.maxloc.amn", U)

[ Info: Written to file: si2.maxloc.amn

Voilà! We have finished the first tutorial! 🎉🎉🎉

This page was generated using Literate.jl.

1. Maximal localization of isolated manifold

Preparation

Tips on initial guess

Model generation

Maximal localization

Save the new gauge