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. Exact exchange does not yet support symmetry-reduced k-point grids, hence we use the full grid already for PBE.
model_PBE = model_DFT(system; pseudopotentials, functionals=PBE())
basis_PBE = PlaneWaveBasis(model_PBE; Ecut=10, kgrid=[3, 3, 3], # Not converged !
use_symmetries_for_kpoint_reduction=false)
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.436107061893 -0.91 5.6 27.8s
2 -8.438609589363 -2.60 -1.78 1.0 7.40s
3 -8.438665009954 -4.26 -2.70 2.4 950ms
4 -8.438668836330 -5.42 -3.54 3.1 462ms
5 -8.438668912654 -7.12 -4.65 3.0 454ms
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,
use_symmetries_for_kpoint_reduction=false)
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.310625453065 -2.16 1.0 9.64s
2 -8.310675928622 -4.30 -2.60 1.0 8.26s
3 -8.310679198856 -5.49 -3.14 11.0 8.92s
4 -8.310679222924 -7.62 -3.82 7.2 8.68s
5 -8.310679223838 -9.04 -4.59 7.4 9.11s
6 -8.310679224112 -9.56 -5.41 9.4 8.83s
7 -8.310679224129 -10.78 -6.20 11.2 9.00s
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.310679224305 -7.36 8.9 25.2s
2 -8.310679224304 + -12.30 -7.88 4.7 27.8s
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.
dis_froz_max = scfres.εF + austrip(2u"eV")
wannier_model = Wannier.Model(scfres;
n_wannier=length(projections),
n_bands=12,
projections,
dis_froz_max)
# 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 # hideb-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: Initial spread
WF center [rx, ry, rz]/Šspread/Ų
1 0.23688 0.23611 0.23703 1.92306
2 0.23710 -0.23709 -0.23676 1.92381
3 -0.23525 0.23514 -0.23557 1.91980
4 -0.23629 -0.23694 0.23710 1.92332
5 1.74045 1.74096 1.74035 1.69101
6 1.74017 0.97482 0.97460 1.69245
7 0.97291 1.74217 0.97312 1.68475
8 0.97421 0.97463 1.74026 1.69150
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.51576
Ω̃ = 4.93394
ΩOD = 4.85682
ΩD = 0.07712
Ω = 14.44970
[ Info: Initial spread (with states freezed)
WF center [rx, ry, rz]/Šspread/Ų
1 0.24231 0.24231 0.24231 1.94248
2 0.24231 -0.24231 -0.24231 1.94248
3 -0.24231 0.24231 -0.24231 1.94248
4 -0.24231 -0.24231 0.24231 1.94248
5 1.72517 1.72517 1.72517 1.75897
6 1.72517 0.98983 0.98983 1.75897
7 0.98983 1.72517 0.98983 1.75897
8 0.98983 0.98983 1.72517 1.75897
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.95815
Ω̃ = 4.84765
ΩOD = 4.77514
ΩD = 0.07251
Ω = 14.80580
┌ 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.480580e+01 6.865670e-01
* time: 0.0012431144714355469
1 1.338353e+01 6.583002e-02
* time: 0.5425190925598145
2 1.301391e+01 4.344046e-02
* time: 0.5485291481018066
3 1.292972e+01 1.919173e-02
* time: 0.5545980930328369
4 1.290488e+01 1.020432e-02
* time: 0.5630450248718262
5 1.289482e+01 5.891568e-03
* time: 0.5740170478820801
6 1.288967e+01 6.661145e-03
* time: 0.5873379707336426
7 1.288663e+01 3.500178e-03
* time: 0.6006171703338623
8 1.288531e+01 2.627752e-03
* time: 0.6139681339263916
9 1.288481e+01 1.531479e-03
* time: 0.6276299953460693
10 1.288470e+01 1.101739e-03
* time: 0.6371080875396729
11 1.288464e+01 5.037163e-04
* time: 0.6508979797363281
12 1.288461e+01 3.595659e-04
* time: 0.6640200614929199
13 1.288460e+01 2.778058e-04
* time: 0.6773350238800049
14 1.288459e+01 2.616216e-04
* time: 0.6905419826507568
15 1.288459e+01 1.970229e-04
* time: 0.7630550861358643
16 1.288458e+01 1.990018e-04
* time: 0.771373987197876
17 1.288458e+01 2.526399e-04
* time: 0.7796871662139893
18 1.288457e+01 2.945123e-04
* time: 0.7879199981689453
19 1.288457e+01 1.322064e-04
* time: 0.7982470989227295
20 1.288457e+01 1.304139e-04
* time: 0.8115551471710205
21 1.288457e+01 1.180300e-04
* time: 0.8247661590576172
[ Info: Final spread
WF center [rx, ry, rz]/Šspread/Ų
1 0.29802 0.29723 0.29623 1.71684
2 0.29462 -0.29602 -0.29489 1.71940
3 -0.30630 0.30606 -0.30592 1.70540
4 -0.29626 -0.29605 0.29713 1.71772
5 1.83119 1.83812 1.83204 1.50765
6 1.83171 0.88411 0.87852 1.50957
7 0.87764 1.83796 0.87884 1.49968
8 0.87726 0.88360 1.83143 1.50830
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.79728
Ω̃ = 3.08729
ΩOD = 2.99089
ΩD = 0.09640
Ω = 12.88457
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.