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.436108520431 -0.91 5.7 27.0s
2 -8.438609933661 -2.60 -1.78 1.0 7.04s
3 -8.438665015907 -4.26 -2.70 2.5 519ms
4 -8.438668835356 -5.42 -3.53 3.3 499ms
5 -8.438668911726 -7.12 -4.51 3.1 998ms
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.310625453359 -2.16 1.0 9.41s
2 -8.310675928876 -4.30 -2.60 1.0 8.28s
3 -8.310679198859 -5.49 -3.14 10.9 9.05s
4 -8.310679223120 -7.62 -3.82 6.6 8.84s
5 -8.310679224276 -8.94 -4.59 5.9 9.22s
6 -8.310679224332 -10.25 -5.41 6.1 8.64s
7 -8.310679224333 -11.76 -6.20 6.9 8.68s
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 25.3s
2 -8.310679224333 -14.75 -7.86 4.7 27.0s
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.23725 0.23412 0.23979 1.92393
2 0.23655 -0.23238 -0.23927 1.92156
3 -0.23695 0.23411 -0.23996 1.92387
4 -0.23597 -0.23209 0.23880 1.92056
5 1.73996 1.74038 1.74004 1.69290
6 1.74085 0.97307 0.97417 1.68832
7 0.97482 1.74049 0.97494 1.69251
8 0.97353 0.97256 1.74145 1.68598
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.51576
Ω̃ = 4.93387
ΩOD = 4.85684
ΩD = 0.07703
Ω = 14.44962
[ Info: Initial spread (with states freezed)
WF center [rx, ry, rz]/Šspread/Ų
1 0.24237 0.23974 0.24483 1.94252
2 0.24237 -0.23973 -0.24483 1.94251
3 -0.24234 0.23975 -0.24484 1.94245
4 -0.24234 -0.23974 0.24483 1.94241
5 1.72505 1.72389 1.72643 1.75913
6 1.72506 0.99113 0.98854 1.75909
7 0.98994 1.72397 0.98847 1.75889
8 0.98989 0.99101 1.72657 1.75874
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.95815
Ω̃ = 4.84758
ΩOD = 4.77514
ΩD = 0.07244
Ω = 14.80573
┌ 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.480573e+01 6.791598e-01
* time: 0.0012700557708740234
1 1.338743e+01 6.545652e-02
* time: 0.5276570320129395
2 1.301998e+01 4.987689e-02
* time: 0.5336880683898926
3 1.292872e+01 2.502016e-02
* time: 0.5396130084991455
4 1.290293e+01 1.032572e-02
* time: 0.5479881763458252
5 1.289247e+01 6.191572e-03
* time: 0.5594460964202881
6 1.288702e+01 6.631201e-03
* time: 0.5727050304412842
7 1.288385e+01 4.510400e-03
* time: 0.5863091945648193
8 1.288248e+01 2.832606e-03
* time: 0.5999100208282471
9 1.288197e+01 1.426563e-03
* time: 0.6134440898895264
10 1.288184e+01 9.677777e-04
* time: 0.6229250431060791
11 1.288177e+01 6.611875e-04
* time: 0.6365330219268799
12 1.288174e+01 4.661429e-04
* time: 0.6502561569213867
13 1.288173e+01 2.964086e-04
* time: 0.7277271747589111
14 1.288172e+01 2.427587e-04
* time: 0.7362301349639893
15 1.288172e+01 1.219262e-04
* time: 0.7443511486053467
16 1.288171e+01 1.634343e-04
* time: 0.7524471282958984
17 1.288171e+01 2.759306e-04
* time: 0.7607002258300781
18 1.288170e+01 2.424479e-04
* time: 0.771644115447998
19 1.288170e+01 1.595648e-04
* time: 0.7849540710449219
20 1.288170e+01 9.333608e-05
* time: 0.798288106918335
21 1.288170e+01 1.269430e-04
* time: 0.8116250038146973
22 1.288169e+01 1.981001e-04
* time: 0.8249800205230713
23 1.288169e+01 3.918024e-04
* time: 0.838374137878418
24 1.288161e+01 1.619969e-03
* time: 0.8950440883636475
25 1.288146e+01 3.582839e-03
* time: 0.9033210277557373
26 1.288142e+01 3.937142e-03
* time: 0.9091801643371582
27 1.288124e+01 4.888742e-03
* time: 0.9151310920715332
28 1.288103e+01 5.697238e-03
* time: 0.920997142791748
29 1.288042e+01 6.681008e-03
* time: 0.9271330833435059
30 1.287816e+01 7.884847e-03
* time: 0.9404780864715576
31 1.286910e+01 1.308787e-02
* time: 0.9539592266082764
32 1.283823e+01 2.911608e-02
* time: 0.967257022857666
33 1.281851e+01 2.664144e-02
* time: 0.9766850471496582
34 1.273021e+01 4.432211e-02
* time: 0.9860661029815674
35 1.259394e+01 3.086719e-02
* time: 0.995427131652832
36 1.248274e+01 4.289373e-02
* time: 1.0264642238616943
37 1.239443e+01 3.162809e-02
* time: 1.0331032276153564
38 1.230190e+01 3.195026e-02
* time: 1.039161205291748
39 1.224274e+01 3.089826e-02
* time: 1.0475680828094482
40 1.216496e+01 1.790360e-02
* time: 1.0559661388397217
41 1.209926e+01 1.967258e-02
* time: 1.0661170482635498
42 1.204508e+01 1.862355e-02
* time: 1.0795331001281738
43 1.200456e+01 2.228921e-02
* time: 1.0928781032562256
44 1.197099e+01 1.843010e-02
* time: 1.1197152137756348
45 1.195107e+01 1.384270e-02
* time: 1.1281020641326904
46 1.193694e+01 7.960994e-03
* time: 1.1364610195159912
47 1.192764e+01 1.015824e-02
* time: 1.144806146621704
48 1.192137e+01 8.467845e-03
* time: 1.1535820960998535
49 1.191709e+01 4.891461e-03
* time: 1.1670231819152832
50 1.191396e+01 3.778558e-03
* time: 1.187075138092041
51 1.191199e+01 3.017804e-03
* time: 1.1952831745147705
52 1.191070e+01 2.702254e-03
* time: 1.2035880088806152
53 1.190987e+01 2.542900e-03
* time: 1.2117881774902344
54 1.190933e+01 2.280740e-03
* time: 1.2211930751800537
55 1.190886e+01 2.162483e-03
* time: 1.2435550689697266
56 1.190850e+01 1.665610e-03
* time: 1.2519052028656006
57 1.190826e+01 1.616575e-03
* time: 1.2602131366729736
58 1.190810e+01 8.382220e-04
* time: 1.2684650421142578
59 1.190801e+01 7.277504e-04
* time: 1.2768092155456543
60 1.190796e+01 4.909418e-04
* time: 1.2881500720977783
61 1.190792e+01 3.539038e-04
* time: 1.3097190856933594
62 1.190790e+01 3.435789e-04
* time: 1.3178842067718506
63 1.190789e+01 3.053471e-04
* time: 1.326134204864502
64 1.190789e+01 2.482984e-04
* time: 1.3342721462249756
65 1.190788e+01 1.397895e-04
* time: 1.3424220085144043
66 1.190788e+01 1.540585e-04
* time: 1.3557450771331787
67 1.190788e+01 9.545840e-05
* time: 1.3738582134246826
68 1.190787e+01 9.600810e-05
* time: 1.381937026977539
69 1.190787e+01 6.868254e-05
* time: 1.390071153640747
70 1.190787e+01 4.377842e-05
* time: 1.3981940746307373
[ Info: Final spread
WF center [rx, ry, rz]/Šspread/Ų
1 0.45349 0.45693 -0.45088 1.48752
2 -0.45366 -0.45706 -0.45009 1.48891
3 -0.45325 0.45749 0.45046 1.48849
4 0.45371 -0.45666 0.45070 1.48818
5 1.81096 1.80828 1.81439 1.48900
6 1.81095 0.90692 0.90074 1.48914
7 0.90409 1.80821 0.90051 1.48829
8 0.90402 0.90718 1.81463 1.48835
Sum spread: Ω = ΩI + Ω̃, Ω̃ = ΩOD + ΩD
ΩI = 9.65072
Ω̃ = 2.25715
ΩOD = 2.16540
ΩD = 0.09175
Ω = 11.90787
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.