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.08219356942080991Then 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.0006182370306135084
└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/terms/xc.jl:39
n Energy log10(ΔE) log10(Δρ) Magnet |Magn| Diag Δtime
--- --------------- --------- --------- ------ ------ ---- ------
1 -361.3879241704 0.07 1.335 3.438 6.8 4.69s
2 -363.2380475053 0.27 -0.21 0.014 3.623 3.4 8.47s
3 -363.3508083405 -0.95 -0.58 0.000 3.727 3.2 3.43s
4 -363.3889545821 -1.42 -1.17 0.000 3.717 2.6 2.37s
5 -363.3959523600 -2.16 -1.66 0.000 3.681 2.0 2.21s
6 -363.3973107595 -2.87 -2.04 0.000 3.656 1.5 2.45s
7 -363.3976160642 -3.52 -2.29 0.000 3.647 2.5 2.18s
8 -363.3976924645 -4.12 -2.65 0.000 3.647 1.2 1.88s
9 -363.3977065439 -4.85 -2.95 0.000 3.649 2.1 2.74s
10 -363.3977059445 + -6.22 -2.90 -0.000 3.649 1.6 1.86s
11 -363.3977092560 -5.48 -3.20 -0.000 3.648 1.6 1.94s
12 -363.3977089351 + -6.49 -3.17 0.000 3.648 2.0 2.66s
13 -363.3977091225 -6.73 -3.17 -0.000 3.648 2.0 2.15s
14 -363.3977086815 + -6.36 -3.03 -0.000 3.649 1.0 1.80s
15 -363.3977092332 -6.26 -3.01 -0.000 3.649 1.0 2.34s
16 -363.3977093966 -6.79 -3.10 -0.000 3.649 1.0 1.75s
17 -363.3977096329 -6.63 -3.11 -0.000 3.649 1.0 1.76s
18 -363.3977097086 -7.12 -3.13 0.000 3.649 1.0 2.31s
19 -363.3977097333 -7.61 -3.10 0.000 3.649 1.0 1.74s
20 -363.3977097464 -7.88 -3.19 -0.000 3.649 1.0 1.75s
21 -363.3977098094 -7.20 -3.48 -0.000 3.649 1.0 1.76s
22 -363.3977099805 -6.77 -3.61 0.000 3.648 1.0 2.24s
23 -363.3977100009 -7.69 -3.78 0.000 3.648 1.0 1.67s
24 -363.3977100105 -8.02 -3.93 0.000 3.648 1.1 1.70s
25 -363.3977100160 -8.26 -4.15 0.000 3.648 1.4 1.72s
26 -363.3977100168 -9.06 -4.05 0.000 3.648 1.5 2.27s
27 -363.3977100171 -9.56 -4.09 0.000 3.648 1.0 1.66s
28 -363.3977100175 -9.46 -4.37 0.000 3.648 1.0 1.67s
29 -363.3977100175 -10.33 -4.58 0.000 3.648 1.0 2.22s
30 -363.3977100174 + -10.17 -4.67 0.000 3.648 1.0 1.66s
31 -363.3977100177 -9.63 -4.88 0.000 3.648 1.1 1.68s
32 -363.3977100178 -10.11 -5.19 0.000 3.648 1.0 1.68s
33 -363.3977100178 -10.17 -5.52 0.000 3.648 1.9 2.42s
34 -363.3977100178 -10.71 -5.73 0.000 3.648 2.4 2.06s
35 -363.3977100178 -11.40 -5.73 0.000 3.648 1.1 1.69s
36 -363.3977100178 -11.77 -5.73 0.000 3.648 1.0 2.21s
37 -363.3977100179 -11.81 -5.66 0.000 3.648 1.0 1.65s
38 -363.3977100179 -12.20 -5.66 0.000 3.648 1.0 1.67s
39 -363.3977100179 -12.47 -5.56 0.000 3.648 1.0 2.25s
40 -363.3977100179 -12.25 -5.50 0.000 3.648 1.0 1.67s
41 -363.3977100179 -12.94 -5.47 0.000 3.648 1.0 1.67s
42 -363.3977100179 + -12.64 -5.39 0.000 3.648 1.0 1.67s
43 -363.3977100179 -12.47 -5.47 0.000 3.648 1.0 2.23s
44 -363.3977100179 -12.40 -5.79 0.000 3.648 1.0 1.66s
45 -363.3977100179 -12.77 -6.11 0.000 3.648 1.0 1.67s
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.11667611134672567With 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)