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.082192885579363Then 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.3855800657 0.07 1.334 3.441 6.8 3.81s
2 -363.2388920759 0.27 -0.21 0.014 3.625 3.2 4.95s
3 -363.3509372226 -0.95 -0.58 0.000 3.727 3.2 2.94s
4 -363.3890226901 -1.42 -1.18 0.000 3.717 2.6 2.46s
5 -363.3959840285 -2.16 -1.67 0.000 3.681 2.0 2.27s
6 -363.3973187580 -2.87 -2.04 0.000 3.656 1.5 1.95s
7 -363.3976110674 -3.53 -2.28 0.000 3.648 2.4 3.81s
8 -363.3976915666 -4.09 -2.64 0.000 3.647 1.2 1.98s
9 -363.3977068388 -4.82 -2.98 0.000 3.649 2.1 2.23s
10 -363.3977059365 + -6.04 -2.90 -0.000 3.650 2.1 2.01s
11 -363.3977093110 -5.47 -3.20 0.000 3.648 2.0 2.13s
12 -363.3977081121 + -5.92 -3.08 0.000 3.648 2.1 2.19s
13 -363.3977089829 -6.06 -3.11 -0.000 3.648 2.1 2.24s
14 -363.3977086162 + -6.44 -3.00 -0.000 3.649 1.1 1.81s
15 -363.3977089676 -6.45 -2.93 0.000 3.649 1.0 1.79s
16 -363.3977092547 -6.54 -3.06 -0.000 3.649 1.0 1.81s
17 -363.3977094786 -6.65 -3.07 -0.000 3.649 1.0 1.81s
18 -363.3977096213 -6.85 -3.08 0.000 3.649 1.0 1.79s
19 -363.3977097042 -7.08 -3.09 0.000 3.649 1.0 1.79s
20 -363.3977092213 + -6.32 -3.01 -0.000 3.649 1.0 1.79s
21 -363.3977089449 + -6.56 -2.97 -0.000 3.649 1.0 1.79s
22 -363.3977097080 -6.12 -3.18 0.000 3.649 1.0 1.80s
23 -363.3977098312 -6.91 -3.29 0.000 3.649 1.0 1.79s
24 -363.3977099184 -7.06 -3.35 0.000 3.649 1.0 1.70s
25 -363.3977100101 -7.04 -4.02 0.000 3.648 1.1 1.72s
26 -363.3977100156 -8.26 -4.10 0.000 3.648 2.5 2.02s
27 -363.3977100127 + -8.54 -4.20 0.000 3.648 1.2 1.77s
28 -363.3977100149 -8.67 -4.51 0.000 3.648 1.0 1.70s
29 -363.3977100159 -8.99 -4.57 0.000 3.648 1.0 1.71s
30 -363.3977100153 + -9.23 -4.50 0.000 3.648 1.0 1.69s
31 -363.3977100128 + -8.60 -4.39 0.000 3.648 1.1 1.72s
32 -363.3977100161 -8.48 -4.54 0.000 3.648 1.1 1.72s
33 -363.3977100167 -9.21 -4.51 0.000 3.648 1.0 3.21s
34 -363.3977100178 -8.98 -5.15 0.000 3.648 1.4 1.77s
35 -363.3977100178 -11.00 -5.18 0.000 3.648 2.6 2.07s
36 -363.3977100178 -10.22 -5.51 0.000 3.648 1.0 1.70s
37 -363.3977100178 -10.61 -5.93 0.000 3.648 2.0 2.00s
38 -363.3977100179 -11.43 -6.49 0.000 3.648 2.1 2.10sRun 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.11667613071457816With 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)