Hubbard correction (DFT+U)
In this example, we'll plot the DOS and projected DOS of Nickel Oxide with and without the Hubbard term correction.
using DFTK
using PseudoPotentialData
using Unitful
using UnitfulAtomic
using PlotsDefine the geometry and pseudopotential
a = 7.9 # Nickel Oxide lattice constant in Bohr
lattice = a * [[ 1.0 0.5 0.5];
[ 0.5 1.0 0.5];
[ 0.5 0.5 1.0]]
pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
Ni = ElementPsp(:Ni, pseudopotentials)
O = ElementPsp(:O, pseudopotentials)
atoms = [Ni, O, Ni, O]
positions = [zeros(3), ones(3) / 4, ones(3) / 2, ones(3) * 3 / 4]
magnetic_moments = [2, 0, -1, 0]4-element Vector{Int64}:
2
0
-1
0First, we run an SCF and band computation without the Hubbard term
model = model_DFT(lattice, atoms, positions; temperature=5e-3,
functionals=PBE(), magnetic_moments)
basis = PlaneWaveBasis(model; Ecut=20, kgrid=[2, 2, 2])
scfres = self_consistent_field(basis; tol=1e-6, ρ=guess_density(basis, magnetic_moments))
bands = compute_bands(scfres, MonkhorstPack(4, 4, 4))
lowest_unocc_band = findfirst(ε -> ε-bands.εF > 0, bands.eigenvalues[1])
band_gap = bands.eigenvalues[1][lowest_unocc_band] - bands.eigenvalues[1][lowest_unocc_band-1]0.08219263372874136Then we plot the DOS and the PDOS for the relevant 3D (pseudo)atomic projector
εF = bands.εF
width = 5.0u"eV"
εrange = (εF - austrip(width), εF + austrip(width))
p = plot_dos(bands; εrange, colors=[:red, :red])
plot_pdos(bands; p, iatom=1, label="3D", colors=[:yellow, :orange], εrange)
To perform and Hubbard computation, we have to define the Hubbard manifold and associated constant.
In DFTK there are a few ways to construct the OrbitalManifold. Here, we will apply the Hubbard correction on the 3D orbital of all nickel atoms. To select all nickel atoms, we can:
- Pass the
Nielement directly. - Pass the
:Nisymbol. - Pass the list of atom indices, here
[1, 3].
To select the orbitals, it is recommended to use their label, such as "3D" for PseudoDojo pseudopotentials.
Note that "manifold" is the standard term used in the literature for the set of atomic orbitals used to compute the Hubbard correction, but it is not meant in the mathematical sense.
U = 10u"eV"
# Alternative:
# manifold = OrbitalManifold(:Ni, "3D")
# Alternative:
# manifold = OrbitalManifold([1, 3], "3D")
manifold = OrbitalManifold(Ni, "3D")OrbitalManifold(Ni, "3D")Run SCF with a DFT+U setup, notice the extra_terms keyword argument, setting up the Hubbard +U term.
model = model_DFT(lattice, atoms, positions; extra_terms=[Hubbard(manifold, U)],
functionals=PBE(), temperature=5e-3, magnetic_moments)
basis = PlaneWaveBasis(model; Ecut=20, kgrid=[2, 2, 2])
scfres = self_consistent_field(basis; tol=1e-6, ρ=guess_density(basis, magnetic_moments));┌ Warning: Negative ρcore detected: -0.0006182370306135084
└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/terms/xc.jl:39
n Energy log10(ΔE) log10(Δρ) Magnet |Magn| Diag Δtime
--- --------------- --------- --------- ------ ------ ---- ------
1 -361.3856815546 0.07 1.334 3.442 7.0 3.75s
2 -363.2386087041 0.27 -0.21 0.014 3.625 3.4 4.82s
3 -363.3509789255 -0.95 -0.58 0.000 3.727 3.4 2.94s
4 -363.3890173856 -1.42 -1.18 0.000 3.717 2.6 2.46s
5 -363.3959764395 -2.16 -1.67 0.000 3.681 2.0 2.24s
6 -363.3973163814 -2.87 -2.04 0.000 3.656 1.5 1.93s
7 -363.3976086969 -3.53 -2.28 0.000 3.648 2.2 2.16s
8 -363.3976907286 -4.09 -2.62 0.000 3.647 1.5 1.97s
9 -363.3977068317 -4.79 -2.99 0.000 3.649 2.1 2.20s
10 -363.3977062139 + -6.21 -2.93 -0.000 3.649 2.2 2.06s
11 -363.3977093255 -5.51 -3.19 0.000 3.648 2.0 2.11s
12 -363.3977086775 + -6.19 -3.15 0.000 3.648 2.1 2.12s
13 -363.3977089884 -6.51 -3.12 -0.000 3.648 2.0 2.14s
14 -363.3977085273 + -6.34 -2.99 -0.000 3.648 1.0 1.77s
15 -363.3977090557 -6.28 -2.98 -0.000 3.649 1.0 1.77s
16 -363.3977093809 -6.49 -3.10 -0.000 3.649 1.0 1.77s
17 -363.3977095841 -6.69 -3.13 -0.000 3.649 1.0 1.77s
18 -363.3977097818 -6.70 -3.24 0.000 3.649 1.0 1.77s
19 -363.3977099012 -6.92 -3.38 0.000 3.648 1.0 1.74s
20 -363.3977099422 -7.39 -3.44 0.000 3.648 1.0 1.68s
21 -363.3977099862 -7.36 -3.56 0.000 3.648 1.0 2.90s
22 -363.3977100143 -7.55 -4.19 0.000 3.648 1.2 1.68s
23 -363.3977100138 + -9.30 -4.24 0.000 3.648 2.2 2.02s
24 -363.3977100126 + -8.91 -4.22 0.000 3.648 1.1 1.69s
25 -363.3977099997 + -7.89 -4.05 0.000 3.648 1.8 1.86s
26 -363.3977099952 + -8.35 -4.06 0.000 3.648 1.0 1.68s
27 -363.3977099926 + -8.59 -4.04 0.000 3.648 1.0 1.68s
28 -363.3977099957 -8.51 -4.05 0.000 3.648 1.0 1.67s
29 -363.3977100033 -8.12 -4.00 0.000 3.648 1.0 1.68s
30 -363.3977100054 -8.66 -4.05 0.000 3.648 1.0 1.68s
31 -363.3977100121 -8.18 -4.30 0.000 3.648 1.0 1.67s
32 -363.3977100163 -8.37 -4.61 0.000 3.648 1.2 1.71s
33 -363.3977100169 -9.24 -4.69 0.000 3.648 1.2 1.70s
34 -363.3977100177 -9.11 -4.80 0.000 3.648 1.9 1.86s
35 -363.3977100178 -9.86 -5.54 0.000 3.648 1.2 1.71s
36 -363.3977100178 -10.59 -5.91 0.000 3.648 3.1 2.30s
37 -363.3977100178 -11.33 -6.03 0.000 3.648 2.0 1.82sRun band computation
bands_hub = compute_bands(scfres, MonkhorstPack(4, 4, 4))
lowest_unocc_band = findfirst(ε -> ε-bands_hub.εF > 0, bands_hub.eigenvalues[1])
band_gap = bands_hub.eigenvalues[1][lowest_unocc_band] - bands_hub.eigenvalues[1][lowest_unocc_band-1]0.1166760235072059With the electron localization introduced by the Hubbard term, the band gap has now opened, reflecting the experimental insulating behaviour of Nickel Oxide.
εF = bands_hub.εF
εrange = (εF - austrip(width), εF + austrip(width))
p = plot_dos(bands_hub; p, colors=[:blue, :blue], εrange)
plot_pdos(bands_hub; p, iatom=1, label="3D", colors=[:green, :purple], εrange)