|
| 1 | +# # Hubbard correction (DFT+U) |
| 2 | +# In this example, we'll plot the DOS and projected DOS of Nickel Oxide |
| 3 | +# with and without the Hubbard term correction. |
| 4 | + |
| 5 | +using DFTK |
| 6 | +using PseudoPotentialData |
| 7 | +using Unitful |
| 8 | +using UnitfulAtomic |
| 9 | +using Plots |
| 10 | + |
| 11 | +# Define the geometry and pseudopotential |
| 12 | +a = 7.9 # Nickel Oxide lattice constant in Bohr |
| 13 | +lattice = a * [[ 1.0 0.5 0.5]; |
| 14 | + [ 0.5 1.0 0.5]; |
| 15 | + [ 0.5 0.5 1.0]] |
| 16 | +pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf") |
| 17 | +Ni = ElementPsp(:Ni, pseudopotentials) |
| 18 | +O = ElementPsp(:O, pseudopotentials) |
| 19 | +atoms = [Ni, O, Ni, O] |
| 20 | +positions = [zeros(3), ones(3) / 4, ones(3) / 2, ones(3) * 3 / 4] |
| 21 | +magnetic_moments = [2, 0, -1, 0] |
| 22 | + |
| 23 | +# First, we run an SCF and band computation without the Hubbard term |
| 24 | +model = model_DFT(lattice, atoms, positions; temperature=5e-3, |
| 25 | + functionals=PBE(), magnetic_moments) |
| 26 | +basis = PlaneWaveBasis(model; Ecut=20, kgrid=[2, 2, 2]) |
| 27 | +scfres = self_consistent_field(basis; tol=1e-6, ρ=guess_density(basis, magnetic_moments)) |
| 28 | +bands = compute_bands(scfres, MonkhorstPack(4, 4, 4)) |
| 29 | +lowest_unocc_band = findfirst(ε -> ε-bands.εF > 0, bands.eigenvalues[1]) |
| 30 | +band_gap = bands.eigenvalues[1][lowest_unocc_band] - bands.eigenvalues[1][lowest_unocc_band-1] |
| 31 | + |
| 32 | +# Then we plot the DOS and the PDOS for the relevant 3D (pseudo)atomic projector |
| 33 | +εF = bands.εF |
| 34 | +width = 5.0u"eV" |
| 35 | +εrange = (εF - austrip(width), εF + austrip(width)) |
| 36 | +p = plot_dos(bands; εrange, colors=[:red, :red]) |
| 37 | +plot_pdos(bands; p, iatom=1, label="3D", colors=[:yellow, :orange], εrange) |
| 38 | + |
| 39 | +# To perform and Hubbard computation, we have to define the Hubbard manifold and associated constant. |
| 40 | +# |
| 41 | +# In DFTK there are a few ways to construct the `OrbitalManifold`. |
| 42 | +# Here, we will apply the Hubbard correction on the 3D orbital of all nickel atoms. |
| 43 | +# |
| 44 | +# Note that "manifold" is the standard term used in the literature for the set of atomic orbitals |
| 45 | +# used to compute the Hubbard correction, but it is not meant in the mathematical sense. |
| 46 | +U = 10u"eV" |
| 47 | +manifold = OrbitalManifold(atoms, Ni, "3D") |
| 48 | + |
| 49 | +# Run SCF with a DFT+U setup, notice the `extra_terms` keyword argument, setting up the Hubbard +U term. |
| 50 | +model = model_DFT(lattice, atoms, positions; extra_terms=[Hubbard(manifold, U)], |
| 51 | + functionals=PBE(), temperature=5e-3, magnetic_moments) |
| 52 | +basis = PlaneWaveBasis(model; Ecut=20, kgrid=[2, 2, 2]) |
| 53 | +scfres = self_consistent_field(basis; tol=1e-6, ρ=guess_density(basis, magnetic_moments)) |
| 54 | + |
| 55 | +# Run band computation |
| 56 | +bands_hub = compute_bands(scfres, MonkhorstPack(4, 4, 4)) |
| 57 | +lowest_unocc_band = findfirst(ε -> ε-bands_hub.εF > 0, bands_hub.eigenvalues[1]) |
| 58 | +band_gap = bands_hub.eigenvalues[1][lowest_unocc_band] - bands_hub.eigenvalues[1][lowest_unocc_band-1] |
| 59 | + |
| 60 | +# With the electron localization introduced by the Hubbard term, the band gap has now opened, |
| 61 | +# reflecting the experimental insulating behaviour of Nickel Oxide. |
| 62 | +εF = bands_hub.εF |
| 63 | +εrange = (εF - austrip(width), εF + austrip(width)) |
| 64 | +p = plot_dos(bands_hub; p, colors=[:blue, :blue], εrange) |
| 65 | +plot_pdos(bands_hub; p, iatom=1, label="3D", colors=[:green, :purple], εrange) |
0 commit comments