A gentle introduction to DFT
In [10]:
include("/Users/harrisonlabollita/software/AtomDFT/src/AtomDFT.jl")
WARNING: replacing module AtomDFT.
Out[10]:
Main.AtomDFT
In [20]:
using DFTK, Plots, LinearAlgebra
A gentle introduction to density-functional theory (DFT)¶
Harry LaBolllita
(Botana group, computational condensed matter physics)
Grad2Grad - August 27, 2021
references depending on interest
the problem...¶
Exercise to the reader: solve the many-body Hamiltonian,
$$ \mathcal{H} = \int d{\bf r} \Psi^{\dagger}({\bf r}) \Big [ -\frac{1}{2} \nabla^{2} + V_{\mathrm{ext}} \, \, ({\bf r}) \Big ] \Psi({\bf r}) + \frac{1}{2} \int d{\bf r}d{\bf r}' \Psi^{\dagger}({\bf r})\Psi^{\dagger}({\bf r}')v_{c}\, ({\bf r} - {\bf r}') \Psi({\bf r}')\Psi({\bf r})$$units: $e^{2} = \hbar = m_{e} = 1$ (energies are in Hartree and lengths are in Bohr radii)
Why bother?¶
In principle, we would know everything about the system.
If our scheme is computationally fast, then we can determine bond lengths, lattice structures, band gaps, effects of dopants, etc.
Exercise to the reader: solve the many-body Hamiltonian,
$$ \mathcal{H} = \int d{\bf r} \Psi^{\dagger}({\bf r}) \Big [ -\frac{1}{2} \nabla^{2} + V_{\mathrm{ext}} \, ({\bf r}) \Big ] \Psi({\bf r}) + \frac{1}{2} \int d{\bf r}d{\bf r}' \Psi^{\dagger}({\bf r})\Psi^{\dagger}({\bf r}')v_{c}\, ({\bf r} - {\bf r}') \Psi({\bf r}')\Psi({\bf r})$$units: $e^{2} = \hbar = m_{e} = 1$ (energies are in Hartree and lengths are in Bohr radii)
$$ E = \mathrm{min}_{\Psi} \, \, \langle \Psi | \mathcal{H} | \Psi \rangle $$key idea¶
Any property of the system of interacting particles can be viewed as functional of the ground state density $n({\bf r})$ !
$$ \langle \Psi[n] | \hat{\mathcal{O}} | \Psi[n] \rangle = O[n] $$the density part...¶
Hohenberg and Kohn, PRB 136, B864 (1964)
There exists a mapping from the ground state electron density $n({\bf r})$ to a unique external potential:
$$ n({\bf r}) \rightarrow V_{\mathrm{ext}} \, ({\bf r})$$HK tell us to do the following)
$$ F[n] = \mathrm{min}_{\Psi\rightarrow n} \, \, \, \, \, \langle \Psi | T[n] + E_{\mathrm{int}}\, \, [n] | \Psi \rangle $$$$ E = \mathrm{min}_{n} \, \, \Big ( F[n] + \int d{\bf r} V_{\mathrm{ext}}\,({\bf r}) \, \, n({\bf r}) \Big ),$$HK tell us to do the following
$$ F[n] = \mathrm{min}_{\Psi\rightarrow n} \, \, \, \, \, \langle \Psi | T[n] + E_{\mathrm{int}}\, \, [n] | \Psi \rangle, $$$$ E = \mathrm{min}_{n} \, \, \Big ( F[n] + \int d{\bf r} V_{\mathrm{ext}}\,({\bf r}) \, \, n({\bf r}) \Big )$$This only has set the stage we still do not know the form of the $F[n]$ or the mapping from $n({\bf r}) \rightarrow V_{\mathrm{ext}} \, ({\bf r}) $. Given the complexity of the many-body wave function, we are probably going to need an approximation (or many).
Kohn & Sham save the day...¶
Kohn and Sham, PRB 140, A1133 (1965)
The exact ground state density of the interacting system is equal to that of an auxillary non-interacting system!
Kohn & Sham save the day...¶
Kohn and Sham, PRB 140, A1133 (1965)
The exact ground state density of the interacting system is equal to that of an auxillary non-interacting system
$$ E_{\mathrm{KS}}\, \, [n] = \sum_{i\sigma} \int d{\bf r}\phi_{i\sigma}^{*}({\bf r})\Big ( -\frac{1}{2} \nabla^{2} \Big ) \phi_{i\sigma}\, ({\bf r}) + \int d{\bf r} V_{\mathrm{ext}}\,({\bf r})n({\bf r}) + E_{\mathrm{H}}\, \,[n] + E_{xc}\, \,[n] $$where $$ n({\bf r}) = \sum_{i\sigma \in \mathrm{occ.}} \, \, f_{\mathrm{FD}}\, |\phi_{i\sigma}({\bf r})|^{2},$$ $$ E_{\mathrm{H}}\, \, = \frac{1}{2} \int d{\bf r}d{\bf r}' \frac{n({\bf r})n({\bf r}')}{|{\bf r} - {\bf r}'|} $$ $$ E_{xc} \, \,= \, \, ? $$
we've seen somthing like this work before...¶
Weiss mean-field theory for the Ising model does essentially the same thing.
$$ \mathcal{H} = \sum_{\langle ij \rangle } Js_{i}s_{j} - h \sum_{i}s_{i}, $$$$ \mathcal{H}_{\mathrm{eff}} = -(h + zJm) \, \, s_{o} $$$$ m = \tanh(\beta ( h + z Jm) )$$| mean-field | density-functional | |
|---|---|---|
| observable | $m$ | $n({\bf r})$ |
| auxiallary system | spin in effecive field | electrons in effective potential |
| Weiss field | local mean field | Kohn-Sham potential |
XC functional¶
Local density approximation (LDA), $$ E_{xc} \, [n] = \int d{\bf r} \, n({\bf r}) \varepsilon_{xc} \, (n({\bf r})), $$ where $\varepsilon_{xc} \, (n({\bf r}))$ is the energy density of a uniform electron gas at the density $n({\bf r})$.
Can be "exactly" parameterized with quantum monte carlo calculations to obtain an analytic form.
the algorithm¶
- guess electron density $n({\bf r})$.
- calculate the Kohn-Sham potential: $V_{\mathrm{KS}} = V_{\mathrm{ext}} + V_{\mathrm{H}} + V_{xc}$
- solve Kohn-Sham equations
- calculate new density and energy.
- Some mixing scheme (improves convergence), then continue 1-4 until convergence.
DFT for an atom¶
Let's code this up for a single atom and compare to the NIST database
In [90]:
r, ρ, pot, Econv = Main.AtomDFT.scf(Z=4, verbose=true);
Found bound state at E = -0.500000000071673 Found bound state at E = -0.12498711431249293 Found bound state at E = -0.049918047596319466 Found bound state at E = -0.12499460664543632 Found bound state at E = -0.051611419760514254 Found bound state at E = -0.053967564424788554 adding state: 0.0, -0.500000000071673 with fermi=1.0 adding state: 0.0, -0.12498711431249293 with fermi=1.0 adding state: 0.0, -0.049918047596319466 with fermi=0.0 Total charge = 4.0000000000000036 e⁻ 1 | Etot[Hartree] = -3.5088422613840984 | ΔE = 7.017684522768197 Found bound state at E = -6.422507214625822 Found bound state at E = -0.7931491874720202 Found bound state at E = -0.1285177993313937 Found bound state at E = -0.7129672234943609 Found bound state at E = -0.09998801056223863 Found bound state at E = -0.04621420163443142 adding state: 0.0, -6.422507214625822 with fermi=1.0 adding state: 0.0, -0.7931491874720202 with fermi=1.0 adding state: 0.0, -0.1285177993313937 with fermi=0.0 Total charge = 4.000000000000009 e⁻ 2 | Etot[Hartree] = -22.530195431164444 | ΔE = 38.04270633956069 Found bound state at E = -3.1791463768352255 Found bound state at E = -0.07018365971502674 Found bound state at E = 0.0195979544012009 Found bound state at E = 0.015714127023241084 Found bound state at E = 0.04279362693638569 Found bound state at E = 0.040769676635219396 adding state: 0.0, -3.1791463768352255 with fermi=1.0 adding state: 0.0, -0.07018365971502674 with fermi=1.0 adding state: 0.0, 0.0195979544012009 with fermi=0.0 Total charge = 4.000000000000005 e⁻ 3 | Etot[Hartree] = -12.232587879145456 | ΔE = 20.595215104037976 Found bound state at E = -4.094486593584075 Found bound state at E = -0.2714307498448706 Found bound state at E = -0.016520170251807707 Found bound state at E = -0.13658254924710128 Found bound state at E = 0.004217978583855864 Found bound state at E = 0.020878934599751098 adding state: 0.0, -4.094486593584075 with fermi=1.0 adding state: 0.0, -0.2714307498448706 with fermi=1.0 adding state: 0.0, -0.016520170251807707 with fermi=0.0 Total charge = 3.9999999999999964 e⁻ 4 | Etot[Hartree] = -15.258519418862875 | ΔE = 6.051863079434838 Found bound state at E = -3.7803022181821015 Found bound state at E = -0.1832054277123853 Found bound state at E = 0.007805014087985686 Found bound state at E = -0.057352067736619314 Found bound state at E = 0.023509597471139947 Found bound state at E = 0.03590169791441133 adding state: 0.0, -3.7803022181821015 with fermi=1.0 adding state: 0.0, -0.1832054277123853 with fermi=1.0 adding state: 0.0, 0.007805014087985686 with fermi=0.0 Total charge = 4.000000000000007 e⁻ 5 | Etot[Hartree] = -14.195146107065868 | ΔE = 2.1267466235940127 Found bound state at E = -3.88074835910626 Found bound state at E = -0.20983776150493788 Found bound state at E = 0.0020211905522915913 Found bound state at E = -0.08048784825297102 Found bound state at E = 0.01922940533434236 Found bound state at E = 0.032816583128342215 adding state: 0.0, -3.88074835910626 with fermi=1.0 adding state: 0.0, -0.20983776150493788 with fermi=1.0 adding state: 0.0, 0.0020211905522915913 with fermi=0.0 Total charge = 4.0000000000000115 e⁻ 6 | Etot[Hartree] = -14.53248490531946 | ΔE = 0.6746775965071841 Found bound state at E = -3.84816361943062 Found bound state at E = -0.20095709256083233 Found bound state at E = 0.004123805721710885 Found bound state at E = -0.07267934373951233 Found bound state at E = 0.020799003592565382 Found bound state at E = 0.03398115277199357 adding state: 0.0, -3.84816361943062 with fermi=1.0 adding state: 0.0, -0.20095709256083233 with fermi=1.0 adding state: 0.0, 0.004123805721710885 with fermi=0.0 Total charge = 3.9999999999999987 e⁻ 7 | Etot[Hartree] = -14.422529170373329 | ΔE = 0.2199114698922635 Found bound state at E = -3.8587245980286555 Found bound state at E = -0.20382759799873992 Found bound state at E = 0.003463932976896908 Found bound state at E = -0.07519474548705667 Found bound state at E = 0.020309354091548257 Found bound state at E = 0.03362119641533269 adding state: 0.0, -3.8587245980286555 with fermi=1.0 adding state: 0.0, -0.20382759799873992 with fermi=1.0 adding state: 0.0, 0.003463932976896908 with fermi=0.0 Total charge = 4.000000000000003 e⁻ 8 | Etot[Hartree] = -14.458165065645586 | ΔE = 0.0712717905445146 Found bound state at E = -3.8552958296038597 Found bound state at E = -0.20289231508662167 Found bound state at E = 0.003681215470804248 Found bound state at E = -0.07437407899214885 Found bound state at E = 0.020470865254866726 Found bound state at E = 0.03374033788295501 adding state: 0.0, -3.8552958296038597 with fermi=1.0 adding state: 0.0, -0.20289231508662167 with fermi=1.0 adding state: 0.0, 0.003681215470804248 with fermi=0.0 Total charge = 3.999999999999993 e⁻ 9 | Etot[Hartree] = -14.446587980995272 | ΔE = 0.023154169300628524 Found bound state at E = -3.856408955602817 Found bound state at E = -0.20319594582456163 Found bound state at E = 0.003610896016549356 Found bound state at E = -0.0746404109049332 Found bound state at E = 0.020418626251015752 Found bound state at E = 0.033701840028933165 adding state: 0.0, -3.856408955602817 with fermi=1.0 adding state: 0.0, -0.20319594582456163 with fermi=1.0 adding state: 0.0, 0.003610896016549356 with fermi=0.0 Total charge = 4.0000000000000036 e⁻ 10 | Etot[Hartree] = -14.45034661046008 | ΔE = 0.007517258929617299 Found bound state at E = -3.856047508135193 Found bound state at E = -0.20309730502197354 Found bound state at E = 0.0036337675362031288 Found bound state at E = -0.07455387460644193 Found bound state at E = 0.020435620413880997 Found bound state at E = 0.03371436881241117 adding state: 0.0, -3.856047508135193 with fermi=1.0 adding state: 0.0, -0.20309730502197354 with fermi=1.0 adding state: 0.0, 0.0036337675362031288 with fermi=0.0 Total charge = 4.0000000000000036 e⁻ 11 | Etot[Hartree] = -14.449126019201042 | ΔE = 0.002441182518076346 Found bound state at E = -3.8561648764378647 Found bound state at E = -0.20312933681360526 Found bound state at E = 0.0036263427858223545 Found bound state at E = -0.0745819747840248 Found bound state at E = 0.020430103966305905 Found bound state at E = 0.0337103022561626 adding state: 0.0, -3.8561648764378647 with fermi=1.0 adding state: 0.0, -0.20312933681360526 with fermi=1.0 adding state: 0.0, 0.0036263427858223545 with fermi=0.0 Total charge = 4.0 e⁻ 12 | Etot[Hartree] = -14.449522372854883 | ΔE = 0.0007927073076814395 Found bound state at E = -3.856126763731061 Found bound state at E = -0.20311893448307805 Found bound state at E = 0.003628754292796047 Found bound state at E = -0.07457284908906578 Found bound state at E = 0.020431895708498796 Found bound state at E = 0.03371162313243313 adding state: 0.0, -3.856126763731061 with fermi=1.0 adding state: 0.0, -0.20311893448307805 with fermi=1.0 adding state: 0.0, 0.003628754292796047 with fermi=0.0 Total charge = 4.000000000000003 e⁻ 13 | Etot[Hartree] = -14.449393664131474 | ΔE = 0.00025741744681795353 Found bound state at E = -3.856139140035399 Found bound state at E = -0.2031223124734317 Found bound state at E = 0.0036279712111513865 Found bound state at E = -0.07457581250579919 Found bound state at E = 0.020431313883592282 Found bound state at E = 0.0337111942115598 adding state: 0.0, -3.856139140035399 with fermi=1.0 adding state: 0.0, -0.2031223124734317 with fermi=1.0 adding state: 0.0, 0.0036279712111513865 with fermi=0.0 Total charge = 3.999999999999997 e⁻ 14 | Etot[Hartree] = -14.449435459755557 | ΔE = 8.359124816692542e-5 Found bound state at E = -3.856135121072467 Found bound state at E = -0.20312121552366033 Found bound state at E = 0.0036282255092456744 Found bound state at E = -0.07457485018014508 Found bound state at E = 0.020431502825928673 Found bound state at E = 0.0337113335006076 adding state: 0.0, -3.856135121072467 with fermi=1.0 adding state: 0.0, -0.20312121552366033 with fermi=1.0 adding state: 0.0, 0.0036282255092456744 with fermi=0.0 Total charge = 4.0 e⁻ 15 | Etot[Hartree] = -14.449421887403705 | ΔE = 2.7144703704351514e-5 Found bound state at E = -3.856136426156777 Found bound state at E = -0.20312157174166443 Found bound state at E = 0.003628142928078369 Found bound state at E = -0.07457516268130794 Found bound state at E = 0.02043144146815485 Found bound state at E = 0.03371128826688688 adding state: 0.0, -3.856136426156777 with fermi=1.0 adding state: 0.0, -0.20312157174166443 with fermi=1.0 adding state: 0.0, 0.003628142928078369 with fermi=0.0 Total charge = 4.000000000000008 e⁻ 16 | Etot[Hartree] = -14.449426294776172 | ΔE = 8.814744933971497e-6 Found bound state at E = -3.8561360023580478 Found bound state at E = -0.20312145606902704 Found bound state at E = 0.003628169743505839 Found bound state at E = -0.0745750612048976 Found bound state at E = 0.020431461392209027 Found bound state at E = 0.0337113029549956 adding state: 0.0, -3.8561360023580478 with fermi=1.0 adding state: 0.0, -0.20312145606902704 with fermi=1.0 adding state: 0.0, 0.003628169743505839 with fermi=0.0 Total charge = 4.0 e⁻ 17 | Etot[Hartree] = -14.449424863579825 | ΔE = 2.862392694424898e-6 Found bound state at E = -3.856136139976679 Found bound state at E = -0.20312149362915438 Found bound state at E = 0.003628161036917138 Found bound state at E = -0.07457509415507721 Found bound state at E = 0.020431454923070098 Found bound state at E = 0.03371129818610427 adding state: 0.0, -3.856136139976679 with fermi=1.0 adding state: 0.0, -0.20312149362915438 with fermi=1.0 adding state: 0.0, 0.003628161036917138 with fermi=0.0 Total charge = 4.00000000000001 e⁻ 18 | Etot[Hartree] = -14.449425328338133 | ΔE = 9.295166165657065e-7 Found bound state at E = -3.856136095281011 Found bound state at E = -0.20312148142611072 Found bound state at E = 0.003628163867313044 Found bound state at E = -0.07457508344919334 Found bound state at E = 0.02043145702615556 Found bound state at E = 0.03371129973689814 adding state: 0.0, -3.856136095281011 with fermi=1.0 adding state: 0.0, -0.20312148142611072 with fermi=1.0 adding state: 0.0, 0.003628163867313044 with fermi=0.0 Total charge = 4.000000000000012 e⁻ 19 | Etot[Hartree] = -14.449425177389571 | ΔE = 3.0189712418859926e-7 Found bound state at E = -3.8561361097955085 Found bound state at E = -0.20312148538878344 Found bound state at E = 0.0036281629481351788 Found bound state at E = -0.07457508692570713 Found bound state at E = 0.020431456343170028 Found bound state at E = 0.033711299233265345 adding state: 0.0, -3.8561361097955085 with fermi=1.0 adding state: 0.0, -0.20312148538878344 with fermi=1.0 adding state: 0.0, 0.0036281629481351788 with fermi=0.0 Total charge = 3.999999999999996 e⁻ 20 | Etot[Hartree] = -14.449425226402202 | ΔE = 9.802526079738527e-8
In [91]:
Main.AtomDFT.plot_scf(r, ρ, pot, Econv)
Let's do a calculation¶
Aluminum metal (Al)
face-centered cubic lattice. Electron configuration $[\mathrm{Ne}]3s^{2}3p^{1}$
In [92]:
a = 7.65339
lattice = diagm(fill(a, 3))
terms_LDA = [Kinetic(), # -1/2 Δ
AtomicLocal(), # part of Vₑₓₜ
AtomicNonlocal(), # part of Vₑₓₜ
PspCorrection(), # part of Vₑₓₜ
Hartree(), # Vₕ
Xc(:lda_xc_teter93), # Vxc
Ewald()] # part of Eₙᵤᵪₗₑₐᵣ
#terms_LDA = [Kinetic(), Hartree(), AtomicLocal()]
Al = ElementPsp(:Al, psp=load_psp("hgh/lda/al-q3"))
atoms = [Al => [[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]]
model = Model(lattice; atoms, temperature=1e-3, terms=terms_LDA)
basis = PlaneWaveBasis(model; Ecut=15, kgrid=[2, 2, 2])
scfres = self_consistent_field(basis);
scfres.energies
n Energy Eₙ-Eₙ₋₁ ρout-ρin α Diag --- --------------- --------- -------- ---- ---- 1 -8.355215069899 NaN 1.34e-01 0.80 10.0 2 -8.356980555949 -1.77e-03 2.57e-02 0.80 1.0 3 -8.357343793999 -3.63e-04 2.29e-03 0.80 4.0 4 -8.357342779167 1.01e-06 1.49e-03 0.80 4.0 5 -8.357345481362 -2.70e-06 2.88e-04 0.80 2.0 6 -8.357345574755 -9.34e-08 5.87e-05 0.80 2.0
Out[92]:
Energy breakdown:
Kinetic 3.5494341
AtomicLocal 1.3990691
AtomicNonlocal 1.5589058
PspCorrection -0.8961586
Hartree 0.0168044
Xc -3.2022680
Ewald -10.7831324
total -8.357345574755
In [83]:
heatmap(scfres.ρ[:, :, 1], c=:Blues)
Out[83]:
In [77]:
bands = plot_bandstructure(scfres, kline_density=10);
plot!(bands, xlims=(0, 1.8), ylims =(-15, 10), size=(900,900), ylabel="ε - εF (eV)")
Computing bands along kpath:
Γ -> X -> W -> K -> Γ -> L -> U -> W -> L -> K and U -> X
Diagonalising Hamiltonian kblocks: 100%|████████████████| Time: 0:00:27
Out[77]:
In [78]:
a = 10.679
lattice = diagm(fill(a, 3))
Ge = ElementPsp(:Ge, psp=load_psp("hgh/lda/ge-q4"))
atoms = [Ge => [[0.0, 0.0, 0.0], [0.25, 0.25, 0.25]]]
model = model_LDA(lattice, atoms, temperature=1e-3)
basis = PlaneWaveBasis(model; Ecut=15, kgrid=[2, 2, 2])
scfres = self_consistent_field(basis);
n Free energy Eₙ-Eₙ₋₁ ρout-ρin α Diag --- --------------- --------- -------- ---- ---- 1 -7.788373851151 NaN 1.71e-01 0.80 6.0 2 -7.793070850625 -4.70e-03 3.68e-02 0.80 1.0 3 -7.793465223638 -3.94e-04 9.30e-03 0.80 2.5 4 -7.793537825304 -7.26e-05 2.76e-03 0.80 3.0 5 -7.793544604516 -6.78e-06 2.86e-04 0.80 2.5 6 -7.793544766439 -1.62e-07 5.14e-05 0.80 3.0
In [79]:
plot_bandstructure(scfres, kline_density=10)
Computing bands along kpath:
Γ and P and Z and Q and Γ and F and P -> ₁ and Q -> ₁ and L and Z
Diagonalising Hamiltonian kblocks: 100%|████████████████| Time: 0:00:31
Out[79]:
Breakdown of DFT¶
Sometimes this auxillary picture breaks down
- band gaps
- strong correlation physics (Mott insulators)
My current research: DFT + DMFT (beyond DFT methods)¶
A better way of dealing with correlated physics. Future Grad2Grad talk!