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.08219343412633134Then 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=[1, 1])
plot_pdos(bands; p, iatom=1, label="3D", colors=[3, 4], ε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. It is also possible to set up multiple manifolds with different U values by passing each pair as a separate entry in the Hubbard constructor (i.e. Hubbard(manifold1 => U1, manifold2 => U2, etc.)) or as two vectors (i.e. Hubbard([manifold1, manifold2, etc.], [U1, U2, etc.])).
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));n Energy log10(ΔE) log10(Δρ) Magnet |Magn| Diag Δtime
--- --------------- --------- --------- ------ ------ ---- ------
1 -361.3882357476 0.07 1.335 3.440 6.9 4.11s
2 -362.9619354048 0.20 -0.10 0.223 3.873 2.6 10.7s
3 -363.1924498688 -0.64 -0.20 0.000 3.776 3.2 2.53s
4 -363.2389275504 -1.33 -0.29 0.000 3.782 2.1 2.05s
5 -363.3700495452 -0.88 -0.30 0.000 3.689 4.0 3.52s
6 -363.3864344127 -1.79 -0.48 -0.000 3.658 2.0 2.05s
7 -363.3967675869 -1.99 -1.13 -0.000 3.676 2.8 2.28s
8 -363.3937249642 + -2.52 -0.92 0.000 3.677 2.0 2.62s
9 -363.3967740008 -2.52 -1.08 0.000 3.656 1.1 1.71s
10 -363.3975246406 -3.12 -1.38 0.000 3.645 1.5 1.85s
11 -363.3976064945 -4.09 -1.47 0.000 3.643 1.0 1.72s
12 -363.3976278760 -4.67 -1.49 0.000 3.643 1.0 2.22s
13 -363.3976772029 -4.31 -1.63 0.000 3.641 1.0 1.65s
14 -363.3974959489 + -3.74 -1.88 -0.000 3.635 1.0 1.65s
15 -363.3976578476 -3.79 -2.21 0.000 3.651 2.0 2.52s
16 -363.3977015430 -4.36 -2.63 -0.000 3.650 1.0 1.65s
17 -363.3976984630 + -5.51 -2.58 -0.000 3.651 1.0 1.63s
18 -363.3976801676 + -4.74 -2.41 -0.000 3.653 1.4 2.27s
19 -363.3976934812 -4.88 -2.55 0.000 3.652 1.0 1.68s
20 -363.3977011823 -5.11 -2.63 0.000 3.651 1.0 1.65s
21 -363.3977090164 -5.11 -2.62 0.000 3.648 1.6 1.78s
22 -363.3977041378 + -5.31 -2.60 0.000 3.650 1.0 2.23s
23 -363.3976963100 + -5.11 -2.54 0.000 3.651 1.4 1.70s
24 -363.3977014460 -5.29 -2.69 0.000 3.651 1.0 1.66s
25 -363.3977071713 -5.24 -2.97 0.000 3.650 1.0 2.23s
26 -363.3977090817 -5.72 -3.22 0.000 3.649 1.1 1.69s
27 -363.3977095931 -6.29 -3.39 -0.000 3.649 1.9 1.91s
28 -363.3977099985 -6.39 -4.06 -0.000 3.648 2.0 2.02s
29 -363.3977099990 -9.26 -4.02 0.000 3.648 2.0 2.55s
30 -363.3977100129 -7.86 -4.34 0.000 3.648 1.0 1.68s
31 -363.3977100159 -8.52 -4.60 0.000 3.648 1.2 1.69s
32 -363.3977100172 -8.90 -4.99 0.000 3.648 2.0 2.57s
33 -363.3977100174 -9.64 -5.19 0.000 3.648 1.6 1.76s
34 -363.3977100176 -9.78 -5.43 0.000 3.648 2.0 1.95s
35 -363.3977100177 -9.99 -5.38 0.000 3.648 2.0 2.58s
36 -363.3977100177 -10.32 -5.37 0.000 3.648 1.0 1.64s
37 -363.3977100178 -10.56 -5.33 0.000 3.648 1.0 1.66s
38 -363.3977100178 + -12.13 -5.03 0.000 3.648 1.0 1.66s
39 -363.3977100178 -10.23 -5.61 0.000 3.648 2.6 2.58s
40 -363.3977100178 -10.79 -5.28 0.000 3.648 1.9 2.02s
41 -363.3977100178 -11.39 -5.24 0.000 3.648 1.0 1.65s
42 -363.3977100179 -11.42 -5.49 0.000 3.648 1.0 2.23s
43 -363.3977100179 -11.83 -5.62 0.000 3.648 1.0 1.68s
44 -363.3977100179 -12.07 -5.72 0.000 3.648 1.0 1.64s
45 -363.3977100179 -12.13 -5.93 0.000 3.648 1.0 1.66s
46 -363.3977100179 -12.40 -6.22 0.000 3.648 1.1 2.24s
Run 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.1166761312083105With 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=[2, 2], εrange)
plot_pdos(bands_hub; p, iatom=1, label="3D", colors=[3, 4], εrange)