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.08219328353577687Then 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.3827031880 0.07 1.334 3.444 6.6 11.3s
2 -362.9704464296 0.20 -0.10 0.224 3.874 2.5 5.51s
3 -363.1933410759 -0.65 -0.20 0.000 3.777 3.1 2.11s
4 -363.2411665048 -1.32 -0.29 0.000 3.781 2.1 2.07s
5 -363.3693277841 -0.89 -0.31 0.000 3.684 4.0 2.36s
6 -363.3861294525 -1.77 -0.48 -0.000 3.653 2.0 1.71s
7 -363.3969226735 -1.97 -1.11 -0.000 3.674 2.8 2.27s
8 -363.3934351424 + -2.46 -0.93 0.000 3.678 2.0 1.71s
9 -363.3967335560 -2.48 -1.08 0.000 3.657 1.0 1.40s
10 -363.3974640925 -3.14 -1.34 0.000 3.647 1.0 1.76s
11 -363.3975730126 -3.96 -1.43 0.000 3.644 1.0 1.46s
12 -363.3975914106 -4.74 -1.42 0.000 3.642 1.0 1.42s
13 -363.3976968665 -3.98 -2.06 0.000 3.643 1.0 1.77s
14 -363.3977050514 -5.09 -2.61 -0.000 3.645 1.0 1.46s
15 -363.3977029832 + -5.68 -2.54 -0.000 3.646 1.6 1.95s
16 -363.3977040712 -5.96 -2.60 -0.000 3.646 1.0 1.48s
17 -363.3977058897 -5.74 -2.74 0.000 3.646 1.0 1.41s
18 -363.3977087042 -5.55 -3.03 0.000 3.647 1.0 1.43s
19 -363.3977097086 -6.00 -3.27 -0.000 3.648 1.6 1.83s
20 -363.3977099541 -6.61 -3.73 0.000 3.648 2.0 1.72s
21 -363.3977098204 + -6.87 -3.46 -0.000 3.649 2.0 1.68s
22 -363.3977099277 -6.97 -3.72 -0.000 3.649 1.0 1.77s
23 -363.3977100038 -7.12 -4.06 0.000 3.648 2.0 1.61s
24 -363.3977099988 + -8.30 -3.97 0.000 3.649 1.0 1.41s
25 -363.3977099982 + -9.20 -3.94 0.000 3.649 1.0 1.76s
26 -363.3977100115 -7.87 -4.08 0.000 3.648 1.0 1.41s
27 -363.3977100161 -8.34 -3.95 0.000 3.648 1.4 1.44s
28 -363.3977100159 + -9.62 -3.92 0.000 3.648 1.0 1.76s
29 -363.3977100162 -9.48 -3.91 0.000 3.648 1.0 1.40s
30 -363.3977100167 -9.30 -3.91 0.000 3.648 1.0 1.41s
31 -363.3977100166 + -10.03 -3.88 0.000 3.648 1.0 1.77s
32 -363.3977100169 -9.63 -3.92 0.000 3.648 1.0 1.40s
33 -363.3977100168 + -10.41 -3.96 0.000 3.648 1.0 1.41s
34 -363.3977100168 + -10.27 -4.01 0.000 3.648 1.0 1.42s
35 -363.3977100166 + -9.66 -4.02 0.000 3.648 1.0 1.77s
36 -363.3977100174 -9.09 -4.35 0.000 3.648 1.0 1.40s
37 -363.3977100177 -9.49 -4.68 0.000 3.648 1.0 1.41s
38 -363.3977100175 + -9.76 -4.88 0.000 3.648 1.1 1.77s
39 -363.3977100177 -9.75 -5.10 0.000 3.648 1.0 1.40s
40 -363.3977100178 -9.91 -5.38 0.000 3.648 2.0 1.59s
41 -363.3977100178 -10.83 -5.68 0.000 3.648 1.0 1.76s
42 -363.3977100178 -10.92 -6.09 0.000 3.648 2.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.11667606559950411With 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)