Hybrid DFT band structures
In this example we discuss how to compute hybrid DFT bands structures using DFTK.
Hybrid-DFT in DFTK is still in its early stages. MPI parallelisation does not work and all treatment of $k$-point symmetry is disabled when using exact exchange terms. We appreciate any issues, bug reports or PRs.
using AtomsBuilder
using DFTK
using PseudoPotentialData
using Wannier
using Brillouin
using CairoMakie
using Unitful
using UnitfulAtomic
pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_5.stringent.upf")
system = bulk(:Si)FlexibleSystem(Si₂, periodicity = TTT):
cell_vectors : [ 0 2.715 2.715;
2.715 0 2.715;
2.715 2.715 0]u"Å"
Atom(Si, [ 0, 0, 0]u"Å")
Atom(Si, [ 1.3575, 1.3575, 1.3575]u"Å")
SCF calculation
We first perform a PBE calculation to get a good starting point. Since we will later Wannierise we make sure to converge 16 orbitals throughout the computational procedure.
model_PBE = model_DFT(system; pseudopotentials, functionals=PBE(), symmetries=false)
basis_PBE = PlaneWaveBasis(model_PBE; Ecut=10, kgrid=[3, 3, 3]) # Not converged !
nbandsalg = AdaptiveBands(model_PBE; n_bands_converge=16)
scfres_PBE = self_consistent_field(basis_PBE; tol=1e-4, nbandsalg)
nothing # hiden Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.436108656017 -0.91 5.6 25.6s
2 -8.438609767295 -2.60 -1.78 1.0 7.19s
3 -8.438665056985 -4.26 -2.70 2.4 377ms
4 -8.438668836241 -5.42 -3.53 3.3 685ms
5 -8.438668912573 -7.12 -4.63 3.1 1.09s
Next we converge an HSE calculation with the default number of bands, again in two steps. First using a more aggressive ACE compression techniques that is more memory efficient and slightly faster (by ignoring all unoccupied orbitals), but that gives the wrong unoccupied bands.
model = model_DFT(system; pseudopotentials, functionals=HSE())
basis = PlaneWaveBasis(model; basis_PBE.Ecut, basis_PBE.kgrid)
exxalg = AceExx(; sketch_with_extra_orbitals=false)
scfres = self_consistent_field(basis; tol=1e-6, scfres_PBE.ρ, scfres_PBE.ψ, scfres_PBE.occupation, exxalg,
scfres_PBE.eigenvalues, solver=ScfDampingSolver(), scfres_PBE.nbandsalg,
diagtolalg=AdaptiveDiagtol(; diagtol_max=1e-2, ratio_ρdiff=0.0005));
nothing # hiden Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.310625453850 -2.16 1.0 9.84s
2 -8.310675928998 -4.30 -2.60 1.0 8.66s
3 -8.310679198861 -5.49 -3.14 10.6 9.45s
4 -8.310679223120 -7.62 -3.82 6.6 9.18s
5 -8.310679224276 -8.94 -4.59 5.9 9.66s
6 -8.310679224332 -10.25 -5.41 6.1 9.06s
7 -8.310679224333 -11.76 -6.20 6.9 9.12s
Finally converge the HSE calculation with laxer ACE compression, that also takes the unoccupied orbitals into account (such that virtual bands are accurate)
exxalg = AceExx(; sketch_with_extra_orbitals=true)
scfres = self_consistent_field(basis; tol=1e-7, scfres.ρ, scfres.ψ, scfres.occupation, exxalg,
scfres.eigenvalues, scfres.solver, scfres.nbandsalg, miniter=2,
diagtolalg=AdaptiveDiagtol(; diagtol_max=1e-6, ratio_ρdiff=0.0005))
nothing # hiden Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.310679224333 -7.37 9.1 26.7s
2 -8.310679224333 -14.75 -7.87 4.7 29.3s
Band structure computation
For the PBE model we can directly use the compute_bands function of DFTK:
kpath = irrfbz_path(model)
bands_PBE = compute_bands(scfres_PBE, kpath; n_bands=12)
# Convert band data to format expected by Brillouin and to eV:
E_PBE = reduce(hcat, bands_PBE.eigenvalues)
E_PBE = ustrip.(auconvert.(u"eV", E_PBE))
E_PBE = E_PBE .- maximum(E_PBE[1:4, :])
nothing # hideFor HSE we proceed via wannierization. For this we setup SP3 projectors for silicon atoms:
projections = [DFTK.wannier_projectors_sp3(model.positions[1]);
DFTK.wannier_projectors_sp3(model.positions[2])]
# Set the frozen window size. Since this is an insulator, about 2eV above the Fermi level
# should be good. Note that Wannier.jl uses eV units for energies:
dis_froz_max = ustrip(auconvert(u"eV", scfres.εF)) + 2.0
wannier_model = Wannier.Model(scfres;
n_wannier=length(projections),
n_bands=12,
projections,
dis_froz_max,
fileprefix="wannier/silicon")
# Run disentanglement, returns a (n_bands, n_wann, n_kpts) array
U = disentangle(wannier_model; max_iter=500);
# InterpModel builds H(R) from model.U, so we store the U giving maximal localisation there
wannier_model.U .= U
nothing # hide[ Info: Reading win file: wannier/silicon.win
num_wann = 8
num_bands = 12
mp_grid = 3 3 3
b-vector shell 1 weight = 0.84022
1 0.38571 -0.38571 -0.38571
2 0.38571 0.38571 -0.38571
3 0.38571 -0.38571 0.38571
4 -0.38571 0.38571 0.38571
5 -0.38571 0.38571 -0.38571
6 0.38571 0.38571 0.38571
7 -0.38571 -0.38571 0.38571
8 -0.38571 -0.38571 -0.38571
Finite difference condition satisfied
[ Info: Reading win file: wannier/silicon.win
num_wann = 8
num_bands = 12
mp_grid = 3 3 3
b-vector shell 1 weight = 0.84022
1 0.38571 -0.38571 -0.38571
2 0.38571 0.38571 -0.38571
3 0.38571 -0.38571 0.38571
4 -0.38571 0.38571 0.38571
5 -0.38571 0.38571 -0.38571
6 0.38571 0.38571 0.38571
7 -0.38571 -0.38571 0.38571
8 -0.38571 -0.38571 -0.38571
Finite difference condition satisfied
[ Info: Reading mmn file: wannier/silicon.mmn
header = Generated by DFTK at 2026-10-30T12:46:10.022
n_bands = 12
n_bvecs = 8
n_kpts = 27
[ Info: Reading wannier/silicon.amn
header = Generated by DFTK at 2026-10-30T12:46:02.822
n_bands = 12
n_wann = 8
n_kpts = 27
[ Info: Reading wannier/silicon.eig
n_bands = 12
n_kpts = 27
[ Info: Initial spread
WF center [rx, ry, rz]/Šspread/Ų
1 0.24051 0.23418 0.23688 1.92412
2 0.23973 -0.23237 -0.23620 1.92151
3 -0.23969 0.23408 -0.23685 1.92366
4 -0.23878 -0.23213 0.23585 1.92058
5 1.73939 1.74010 1.73998 1.69394
6 1.74038 0.97326 0.97408 1.68889
7 0.97485 1.74081 0.97434 1.69147
8 0.97364 0.97234 1.74191 1.68541
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.51576
Ω̃ = 4.93381
ΩOD = 4.85676
ΩD = 0.07706
Ω = 14.44957
[ Info: Initial spread (with states freezed)
WF center [rx, ry, rz]/Šspread/Ų
1 0.24499 0.23975 0.24222 1.94266
2 0.24499 -0.23974 -0.24223 1.94265
3 -0.24486 0.23978 -0.24225 1.94229
4 -0.24486 -0.23977 0.24225 1.94224
5 1.72623 1.72373 1.72497 1.75966
6 1.72625 0.99128 0.99000 1.75962
7 0.98866 1.72418 0.98954 1.75835
8 0.98861 0.99079 1.72549 1.75819
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.95815
Ω̃ = 4.84751
ΩOD = 4.77505
ΩD = 0.07246
Ω = 14.80566
┌ Warning: f_tol is deprecated. Use f_abstol or f_reltol instead. The provided value (1.0e-7) will be used as f_reltol.
└ @ Optim ~/.julia/packages/Optim/gmigl/src/types.jl:120
Iter Function value Gradient norm
0 1.480566e+01 7.178378e-01
* time: 0.0011870861053466797
1 1.338744e+01 6.554645e-02
* time: 0.5312671661376953
2 1.301997e+01 4.520367e-02
* time: 0.5419540405273438
3 1.292871e+01 2.229876e-02
* time: 0.5524911880493164
4 1.290294e+01 1.034070e-02
* time: 0.64894700050354
5 1.289248e+01 6.724239e-03
* time: 0.6579461097717285
6 1.288705e+01 7.222702e-03
* time: 0.6665451526641846
7 1.288389e+01 4.228654e-03
* time: 0.6751711368560791
8 1.288249e+01 2.989752e-03
* time: 0.6836321353912354
9 1.288199e+01 2.142104e-03
* time: 0.6957190036773682
10 1.288181e+01 1.015644e-03
* time: 0.7106971740722656
11 1.288173e+01 7.114473e-04
* time: 0.7255990505218506
12 1.288170e+01 4.142123e-04
* time: 0.7404091358184814
13 1.288169e+01 3.368429e-04
* time: 0.755230188369751
14 1.288168e+01 2.697017e-04
* time: 0.9039270877838135
15 1.288167e+01 2.212987e-04
* time: 0.9128730297088623
16 1.288167e+01 2.764516e-04
* time: 0.9214940071105957
17 1.288166e+01 2.643837e-04
* time: 0.9299380779266357
18 1.288165e+01 2.955413e-04
* time: 0.9383571147918701
19 1.288164e+01 1.846572e-04
* time: 0.9516110420227051
20 1.288164e+01 1.186466e-04
* time: 0.9664139747619629
21 1.288164e+01 1.609347e-04
* time: 0.9811110496520996
22 1.288163e+01 2.849070e-04
* time: 0.9959239959716797
23 1.288157e+01 1.357510e-03
* time: 1.0151221752166748
24 1.288059e+01 8.513220e-03
* time: 1.6519780158996582
25 1.288035e+01 9.563630e-03
* time: 1.658177137374878
26 1.287674e+01 1.202706e-02
* time: 1.664276123046875
27 1.286928e+01 1.794018e-02
* time: 1.67280912399292
28 1.285646e+01 2.141905e-02
* time: 1.678861141204834
29 1.281917e+01 2.151783e-02
* time: 1.689561128616333
30 1.273514e+01 2.593342e-02
* time: 1.7044551372528076
31 1.266499e+01 4.940172e-02
* time: 1.714919090270996
32 1.255946e+01 2.134245e-02
* time: 1.7297801971435547
33 1.248640e+01 1.854899e-02
* time: 1.7448291778564453
34 1.244584e+01 1.870056e-02
* time: 1.7597131729125977
35 1.241912e+01 1.166481e-02
* time: 1.7747001647949219
36 1.239845e+01 1.010310e-02
* time: 1.7896502017974854
37 1.237090e+01 1.653126e-02
* time: 1.8053359985351562
38 1.232540e+01 1.918433e-02
* time: 1.8205959796905518
39 1.225218e+01 1.897774e-02
* time: 1.835815191268921
40 1.216939e+01 1.806427e-02
* time: 1.8508310317993164
41 1.208884e+01 2.331750e-02
* time: 1.8612949848175049
42 1.203662e+01 2.104986e-02
* time: 1.8760809898376465
43 1.200009e+01 2.101621e-02
* time: 1.8909189701080322
44 1.197192e+01 1.060100e-02
* time: 1.9059600830078125
45 1.195630e+01 1.068145e-02
* time: 1.9209530353546143
46 1.194209e+01 6.853367e-03
* time: 1.9358761310577393
47 1.193201e+01 1.054774e-02
* time: 1.9507560729980469
48 1.192506e+01 7.202403e-03
* time: 1.9656531810760498
49 1.191918e+01 5.728085e-03
* time: 1.980696201324463
50 1.191564e+01 4.229356e-03
* time: 1.9957151412963867
51 1.191322e+01 4.065588e-03
* time: 2.0106871128082275
52 1.191168e+01 3.142450e-03
* time: 2.0257561206817627
53 1.191044e+01 2.886366e-03
* time: 2.0406670570373535
54 1.190968e+01 2.533753e-03
* time: 2.0555551052093506
55 1.190908e+01 1.824687e-03
* time: 2.0706350803375244
56 1.190866e+01 1.505069e-03
* time: 2.1825830936431885
57 1.190837e+01 1.081129e-03
* time: 2.1914000511169434
58 1.190816e+01 9.798037e-04
* time: 2.2000529766082764
59 1.190801e+01 7.110444e-04
* time: 2.2085821628570557
60 1.190792e+01 7.273010e-04
* time: 2.2190651893615723
61 1.190787e+01 4.600025e-04
* time: 2.233966112136841
62 1.190785e+01 4.126643e-04
* time: 2.2490031719207764
63 1.190784e+01 3.160198e-04
* time: 2.263942003250122
64 1.190783e+01 1.850761e-04
* time: 2.2789270877838135
65 1.190782e+01 1.463546e-04
* time: 2.2941479682922363
66 1.190782e+01 1.336355e-04
* time: 2.3092870712280273
67 1.190782e+01 1.080243e-04
* time: 2.3244130611419678
68 1.190782e+01 9.382751e-05
* time: 2.3395121097564697
69 1.190782e+01 6.433879e-05
* time: 2.3544411659240723
[ Info: Final spread
WF center [rx, ry, rz]/Šspread/Ų
1 -0.45014 -0.45680 -0.45375 1.49031
2 0.45040 -0.45678 0.45393 1.48664
3 0.45038 0.45703 -0.45430 1.48597
4 -0.44992 0.45745 0.45342 1.48999
5 1.81430 1.80825 1.81074 1.49062
6 1.81386 0.90711 0.90361 1.49100
7 0.90082 1.80848 0.90365 1.48674
8 0.90084 0.90693 1.81130 1.48655
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.65076
Ω̃ = 2.25706
ΩOD = 2.16531
ΩD = 0.09175
Ω = 11.90782
Based on the Wannier functions compute the bands via Wannier interpolation
interp_model = Wannier.InterpModel(wannier_model; kpath)
kinter = Brillouin.interpolate(kpath, density=40) # same density as in DFTK by default
data = Wannier.interpolate(interp_model, kinter)
E_HSE = data .- maximum(data[1:4, :])
nothing # hidePlot the bands
Finally we plot the bands using the Makie interface of Brillouin:
fig = Figure(size=(700, 500))
CairoMakie.plot(fig[1, 1], kinter, permutedims(E_HSE), permutedims(E_PBE);
label=["HSE", "PBE"],
ylabel="E - E_VBM (eV)",
ylims=(-13.0, 10.0),
linewidth=2,
linestyle=[:solid, :dash],
color=[:royalblue, :black])
fig
Further reading
For more details on Wannierisation, see also Wannierization using Wannier.jl or Wannier90.
(Delete temporary files when done.)
rm("wannier", recursive=true)