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.08219377403973471Then 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));┌ Warning: Negative ρcore detected: -0.0006182370306135009
└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/terms/xc.jl:39
n Energy log10(ΔE) log10(Δρ) Magnet |Magn| Diag Δtime
--- --------------- --------- --------- ------ ------ ---- ------
1 -361.3875035617 0.07 1.335 3.439 6.9 4.10s
2 -363.2387743959 0.27 -0.21 0.014 3.624 3.0 10.4s
3 -363.3504473069 -0.95 -0.58 0.000 3.727 3.2 3.45s
4 -363.3889407564 -1.41 -1.17 0.000 3.717 2.8 2.40s
5 -363.3959570778 -2.15 -1.66 0.000 3.681 2.0 2.14s
6 -363.3973141492 -2.87 -2.04 0.000 3.656 1.8 1.90s
7 -363.3976088385 -3.53 -2.28 0.000 3.647 2.1 2.69s
8 -363.3976885783 -4.10 -2.60 0.000 3.647 1.6 1.89s
9 -363.3977053570 -4.78 -2.87 0.000 3.649 2.1 2.09s
10 -363.3977059252 -6.25 -2.88 -0.000 3.649 1.2 2.38s
11 -363.3977075047 -5.80 -2.95 -0.000 3.649 1.0 1.77s
12 -363.3977067301 + -6.11 -2.76 -0.000 3.649 1.0 1.71s
13 -363.3977097291 -5.52 -3.17 0.000 3.648 1.4 1.78s
14 -363.3977098018 -7.14 -3.12 0.000 3.648 1.0 2.30s
15 -363.3977098846 -7.08 -3.11 0.000 3.648 1.0 1.70s
16 -363.3977098722 + -7.91 -3.06 0.000 3.648 1.0 1.68s
17 -363.3977099409 -7.16 -3.22 0.000 3.648 1.0 1.70s
18 -363.3977099879 -7.33 -3.33 0.000 3.648 1.0 2.29s
19 -363.3977100021 -7.85 -3.39 0.000 3.648 1.0 1.63s
20 -363.3977100086 -8.18 -3.53 0.000 3.648 1.0 1.59s
21 -363.3977100113 -8.58 -3.59 0.000 3.648 1.0 1.59s
22 -363.3977099999 + -7.95 -3.62 0.000 3.648 1.0 2.20s
23 -363.3977099773 + -7.64 -3.74 0.000 3.648 1.0 1.59s
24 -363.3977099598 + -7.76 -3.74 -0.000 3.648 2.0 1.91s
25 -363.3977099724 -7.90 -3.77 -0.000 3.649 1.0 1.59s
26 -363.3977100029 -7.52 -4.03 0.000 3.648 1.6 1.83s
27 -363.3977100161 -7.88 -4.75 0.000 3.648 1.5 2.23s
28 -363.3977100162 -9.94 -4.56 -0.000 3.648 3.2 2.33s
29 -363.3977100169 -9.13 -4.74 0.000 3.648 1.0 1.59s
30 -363.3977100177 -9.12 -5.20 0.000 3.648 1.9 1.82s
31 -363.3977100178 -9.94 -5.66 0.000 3.648 2.2 2.59s
32 -363.3977100178 + -10.62 -5.33 0.000 3.648 2.5 2.02s
33 -363.3977100178 + -11.34 -5.24 0.000 3.648 1.0 1.59s
34 -363.3977100178 -10.28 -5.19 0.000 3.648 2.0 1.92s
35 -363.3977100178 -11.18 -5.18 0.000 3.648 1.0 2.21s
36 -363.3977100178 -11.69 -5.15 0.000 3.648 1.0 1.58s
37 -363.3977100178 -12.29 -5.10 0.000 3.648 1.0 1.58s
38 -363.3977100178 -12.94 -5.13 0.000 3.648 1.0 1.60s
39 -363.3977100178 -12.13 -5.18 0.000 3.648 1.0 1.66s
40 -363.3977100178 + -12.77 -5.18 0.000 3.648 1.0 2.18s
41 -363.3977100178 -13.25 -5.11 0.000 3.648 1.0 1.59s
42 -363.3977100178 -12.25 -5.17 0.000 3.648 1.0 1.59s
43 -363.3977100178 -12.77 -5.15 0.000 3.648 1.0 1.59s
44 -363.3977100178 + -10.80 -4.78 0.000 3.648 2.0 2.62s
45 -363.3977100178 + -10.84 -4.74 0.000 3.648 1.0 1.58s
46 -363.3977100179 -10.47 -5.58 0.000 3.648 2.0 2.20s
47 -363.3977100179 -11.94 -5.59 0.000 3.648 1.0 1.59s
48 -363.3977100179 + -12.94 -5.55 0.000 3.648 1.0 2.23s
49 -363.3977100179 + -12.20 -5.38 0.000 3.648 1.0 1.58s
50 -363.3977100179 + -12.34 -5.32 0.000 3.648 1.0 1.58s
51 -363.3977100179 -12.40 -6.10 0.000 3.648 1.0 1.58s
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.11667623215712192With 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)