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.0821933794872553Then 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.3871485422 0.07 1.335 3.440 7.0 4.37s
2 -362.9625131860 0.20 -0.10 0.223 3.873 2.6 10.6s
3 -363.1926582582 -0.64 -0.20 0.000 3.776 3.2 2.99s
4 -363.2395466832 -1.33 -0.29 0.000 3.781 2.2 2.15s
5 -363.3696928849 -0.89 -0.30 0.000 3.689 4.0 3.03s
6 -363.3864120881 -1.78 -0.48 -0.000 3.657 2.0 2.42s
7 -363.3967754722 -1.98 -1.13 -0.000 3.676 2.8 2.35s
8 -363.3935523805 + -2.49 -0.90 0.000 3.677 2.0 2.14s
9 -363.3967723585 -2.49 -1.08 0.000 3.656 1.0 2.01s
10 -363.3975256344 -3.12 -1.39 0.000 3.645 1.5 1.85s
11 -363.3976075591 -4.09 -1.47 0.000 3.643 1.0 2.06s
12 -363.3976312898 -4.62 -1.50 0.000 3.643 1.0 1.72s
13 -363.3976717465 -4.39 -1.59 0.000 3.641 1.0 1.69s
14 -363.3976726658 -6.04 -2.10 0.000 3.651 1.0 1.70s
15 -363.3976737997 -5.95 -2.36 -0.000 3.653 1.0 2.04s
16 -363.3976857043 -4.92 -2.46 -0.000 3.652 1.0 1.69s
17 -363.3976952770 -5.02 -2.57 -0.000 3.652 1.0 1.70s
18 -363.3977088054 -4.87 -3.15 -0.000 3.650 1.4 2.11s
19 -363.3977098584 -5.98 -3.56 -0.000 3.648 2.6 2.22s
20 -363.3977097871 + -7.15 -3.46 -0.000 3.648 1.6 1.94s
21 -363.3977098029 -7.80 -3.46 -0.000 3.648 1.0 2.08s
22 -363.3977099465 -6.84 -3.77 0.000 3.648 1.0 1.69s
23 -363.3977100053 -7.23 -4.25 0.000 3.648 1.4 1.75s
24 -363.3977100148 -8.02 -4.69 0.000 3.648 2.5 2.50s
25 -363.3977100164 -8.81 -4.74 0.000 3.648 1.8 1.87s
26 -363.3977100170 -9.22 -5.28 0.000 3.648 1.1 1.72s
27 -363.3977100173 -9.45 -4.74 0.000 3.648 3.0 2.64s
28 -363.3977100176 -9.67 -4.96 0.000 3.648 1.0 1.71s
29 -363.3977100177 -9.92 -5.96 0.000 3.648 1.6 1.89s
30 -363.3977100177 -10.14 -5.45 0.000 3.648 3.4 2.84s
31 -363.3977100178 -10.38 -5.80 0.000 3.648 2.8 2.38s
32 -363.3977100178 -10.62 -5.44 0.000 3.648 2.0 2.32s
33 -363.3977100178 -10.80 -5.89 0.000 3.648 1.4 2.17s
34 -363.3977100178 -11.05 -5.51 0.000 3.648 1.6 1.87s
35 -363.3977100178 -11.27 -5.43 0.000 3.648 1.4 1.81s
36 -363.3977100178 -11.56 -5.31 0.000 3.648 1.1 1.78s
37 -363.3977100178 -11.81 -5.22 0.000 3.648 1.0 1.99s
38 -363.3977100179 -11.92 -5.26 0.000 3.648 1.0 1.69s
39 -363.3977100179 -12.13 -5.26 0.000 3.648 1.0 1.76s
40 -363.3977100179 -12.40 -5.24 0.000 3.648 1.0 2.02s
41 -363.3977100179 + -12.17 -5.16 0.000 3.648 1.0 1.72s
42 -363.3977100179 -11.99 -5.28 0.000 3.648 1.0 1.71s
43 -363.3977100179 -12.64 -5.28 0.000 3.648 1.0 2.09s
44 -363.3977100179 -11.94 -5.72 0.000 3.648 1.0 1.70s
45 -363.3977100179 -12.34 -6.17 0.000 3.648 1.0 1.71s
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.11667610107634002With 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)