Introduction to DFT

Density Functional Theory (DFT)

DFT is a theoretical framework based on the Hohenberg and Kohn theorems, that established a one-to-one correspondence between the ground state electron density, \(n_0(\mathbf{r})\), the ground state wavefunction, \(\Psi_0(\mathbf{r})\), and the external potential, \(v(\mathbf{r})\). This means that the ground state properties of a system can be fully determined by the density.

\[\psi_0 \leftrightarrow n_0 \leftrightarrow v(\mathbf{r})\]

The ground state total energy is then given as a functional of the electron density, \(E[n]\), which is given by:

\[E[n] = T[n] + E_{ee}[n] + E_{eN}[n] + E_{NN}[n]\]

where \(T[n]\) is the kinetic energy, \(E_{ee}[n]\) is the electron-electron interaction, \(E_{eN}[n]\) the electron-nuclear interaction, and \(E_{NN}[n]\) is the nuclear-nuclear interaction.

To determine the density of the system and evaluate the energy functional, DFT explores the mapping between the real system and fictitious systems shown in Figure 1.

_images/maps.svg

Figure1. Mapping between the real system and fictitious systems in DFT.

Kohn-Sham System

In the non-interacting fermionic system, known as the single-particle or Kohn-Sham system, the total energy functional is given by:

\[E[n] = T_s[n] + E_H[n] + E_{xc}[n] + \int v_{eN}(\mathbf{r}) n(\mathbf{r}) d\mathbf{r} + E_{NN}[n]\]

where the density \(n(r) = \sum_i |\phi_i(r)|^2\). And the single-particle kinetic energy is

\[T_{s}[n] \equiv T_s[\{\phi_i\}]= -\frac{1}{2}\sum_i n_i \langle \phi_i | \nabla^2 | \phi_i\rangle = -\frac{1}{2}\sum_i n_i \int \phi_i^*(r) \nabla^2 \phi_i(r) dr\]

The electron-electron repulsion and the total kinetic energy are related to \(T_s\) and \(E_{xc}\) as follows:

\[E_H[n]=\frac{1}{2}\int \frac{n(r)n(r')}{|r-r'|}drdr'\]
\[E_{xc} = T[n] + E_{ee}[n] = T_s[n] + E_H[n] + E_{xc}[n] \to \textbf{Approximated}\]

For the Kohn-Sham system, the Lagrangian is defined as follows:

\[\mathcal{L}_{KS}[\{\phi_i\}] = E[\{\phi_i\}] - \sum_{ij} \varepsilon_{ij}\left(\langle \phi_j|\phi_i \rangle - \delta_{ij}\right)\]

to find the ground state KS orbitals and ground state density, the KS lagrangian is minimized \(\frac{\delta \mathcal{L}_{KS}[\{\phi_i\}]}{\delta \langle \phi_j|}=0\), or just \(\frac{\delta \mathcal{L}_{KS}[\{\phi_i\}]}{\delta \phi_j^*(r)}=0\), which yields to the KS-equation for the non-interacting fermionic system:

Important

\[\left(-\frac{1}{2} \nabla^2 + v_s(\mathbf{r})\right) \phi_i(\mathbf{r}) = \epsilon_i \phi_i(\mathbf{r})\]

where \(\psi_i(\mathbf{r})\) are the KS orbitals, and \(v_s(\mathbf{r})\) is the Kohn-Sham potential given by:

\[v_s[n](r) = \frac{\delta E_{H}[n]}{\delta n(r)} + \frac{\delta E_{xc}[n]}{\delta n(r)} + v_{eN}(r)\]

Orbital-Free DFT

In the non-interacting bosonic system, used in Orbital-Free DFT, the total energy functionals is given by:

\[E[n] = T_{B}[n] + T_P[n] + E_{H}[n] + E_{xc}[n] + \int v_{eN}[n](r) n(r) dr + E_{NN}\]

where \(T_P[n] = T_{s}[n] - T_{vW}[n]\), \(T_{s}[n]\) is the non-interacting kinetic energy functional (KEDF, which is an \(\textbf{approximated}\) functional), and \(T_{B}[n]\) is the bosonic kinetic energy functional (von Weizsäcker kinetic energy functional), given by:

\[T_{B}[n] = T_{vW}[n] = -\frac{1}{2}\int \phi^*(r) \nabla^2 \phi(r) dr\]

where: \(\phi(r)=\sqrt{n(r)}\). The appropiate Lagrangian for OF-DFT is given by:

\[\mathcal{L}_{OF}[n] = E[n] - \mu \left( \int n(\mathbf{r}) d\mathbf{r} - N \right)\]

where \(\mu\) is the Lagrange multiplier, \(N\) is the number of valence electrons in the system. Both \(n_0(\mathbf{r})\) and \(\mu\) are determined during the minimization.

To find the ground state density, the Lagrangian is minimized \(\frac{\delta \mathcal{L}_{OF}[n]}{\delta \langle \phi|}=0\) or \(\frac{\delta \mathcal{L}_{OF}[n]}{\delta \phi^*(r)}=0\), which yields to the KS-equation for the non-interacting bosonic system:

Important

\[-\frac{1}{2}\nabla^2 \phi(r) + v_B[n](r)\phi(r) = \mu\phi(r)\]

where \(\phi(r)\) is the bosonic orbitals, and \(v_B[n](r)\) is the bosonic potential given by:

\[v_B[n](r) = \underbrace{\frac{\delta T_{s}[n]}{\delta n(r)} - \frac{\delta T_{vW}[n]}{\delta n(r)}}_{\frac{\delta T_P[n]}{\delta n(r)}} + \frac{\delta E_{H}[n]}{\delta n(r)} + \frac{\delta E_{xc}[n]}{\delta n(r)} + v_{eN}(r)\]

Solvers for OF-DFT and KS-DFT

In KS-DFT, the ground state electron density, \(n_0(\mathbf{r})\), is obtained from a self-consistent field (SCF) iteration, as shown in the table below.

In OF-DFT, the ground state electron density, \(n_0(\mathbf{r})\), is obtained from the direct minimization of the ground state energy density functional, \(E[n]\), as shown in the table below. The ground state energy is, \(E_0 = E[n_0]\).

Method

OF-DFT

KS-DFT

Direct Minimization

\(n_0(\mathbf{r})=\arg\underset{n}{\min}\big\{ \mathcal{L}_{OF}[n]\big\}\)

\(\{\phi_i^0\} = \arg\underset{\{\phi_i\}}{\min}\big\{ \mathcal{L}_{KS}[\{\phi_i\}]\big\}\)

SCF

\(-\frac{1}{2}\nabla^2 \sqrt{n(r)} + v_B(r)\sqrt{n(r)} = \mu\sqrt{n(r)}\)

\(-\frac{1}{2}\nabla^2 \phi_i(r) + v_s(r)\phi_i(r) = \varepsilon_i\phi_i(r)\)

Note

In DFTpy, \(T_s[n]\) and \(E_{xc}[n]\) are pure functionals of the density. Check out the tutorials for a list of available KEDF and EXC functionals,

Note

DFTpy solves the ground state problem with the so-called direct energy minimization. Other (faster) methods are available, such as OESCF, which is implemented in eDFTpy. OESCF may be implemented in DFTpy upon request.

Time-Dependent Density-Functional Theory (TDDFT)

DFT can also describe a system out of equilibrium by propagating it in real time or in frequency space finding the roots of the frequency dependent polarizability (Casida). In the time regime the Runge and Gross theorem is an analog of the Hohenberg and Khon theorem. This theorem formally defines a one-to-one correspondance between densities and potentials, for any fixed initial many-body state, the it follows that the time-dependent density is a unique functional of the potentials and vice versa. This means that the many-body Hamiltonian \(\hat{H}(t)\) and thus the many-body wave function \(\Psi(t)\) are functionals of \(n(\mathbf{r},t)\) as well

Following the Runge and Gross theorem, we can write time-dependent Schrodinger like equation, namely:

\[\hat{H}(t) \Psi(\mathbf{r},t) = i \frac{d}{dt}\Psi(\mathbf{r},t)\]

In the KS formalism the KS orbitals are propagated with a time-dependent Hamiltonian, given by:

\[\hat{H}(t) = -\frac{1}{2} \nabla^2 + v_s[n](\mathbf{r},t)\]

where the time-dependent KS potential is given by:

\[\]

v_s[n](mathbf r,t) = v_{eN}(mathbf r,t) + v_H[n](mathbf r,t) + v_{xc}[n](mathbf r,t).

In the OF formalism the bosonic wavefunction \(\Psi(\mathbf{r},t)\) is propagated with a time-dependent KS-like Hamiltonian, given by:

\[\hat{H}(t) = -\frac{1}{2} \nabla^2 + v_B[n](\mathbf{r},t).\]

The Bosonic KS-like potential is given by two major contributions

\[v_B(\mathbf{r},t) = v_s[n](\mathbf{r},t) + v_P[n](\mathbf{r},t)\]

where the adiabatic approximation has been invoked fro the xc potential \(v_{xc}[n(t)]\). The Pauli potential is given by adiabatic and nonadiabatic contributions,

\[v_P(\mathbf{r},t)=v_P^{ad}(\mathbf{r},t)+v_P^{nad}[n](\mathbf{r},t).\]

Note

The adiabatic Pauli potential can be specified according to any given KEDF available in DFTpy.

The nonadiabatic contribution is usually neglected in the literature. In DFTpy the JP functional is available,

\[v_P^{nad}(\mathbf{r},t) = -\frac{\pi^3}{12}\left(\frac{6}{k_F^2(\mathbf{r},t)}\mathcal{F}^{-1}\left\{i\mathbf{q}\cdot\mathbf{j}(\mathbf{q},t)\frac{1}{q}\right\}+\frac{1}{k_F^4(\mathbf{r},t)}\mathcal{F}^{-1}\left\{i\mathbf{q}\cdot\mathbf{j}(\mathbf{q},t)q\right\}\right)\]

where \(\mathbf{j}\) and \(\mathbf{q}\) are the electronic current density and the reciprocal space vector, respectively. The current density is determined by the standard equation \(\mathbf{j}(\mathbf{r})=\frac{1}{2i}\left[\Psi^*(\mathbf{r})\nabla\Psi(\mathbf{r})-\Psi(\mathbf{r})\nabla\Psi^*(\mathbf{r})\right]\). \(\mathcal{F}\) stands for Fourier transform and \(k_F(\mathbf{r},t)=[3\pi^2 n(\mathbf{r},t)]^{1/3}\) is the Fermi wavevector function of the local electron density.

Warning

The JP potential is numerically challenging. Refer to the original JP publication for details.

Note

Optical spectra and nonlinear electronic processes can be modelled by DFTpy. See tutorials for additional information. Ehrenfest dynamics is not yet available.

Linear response theory

Linear-response TDDFT refers to the determination of excitation energies and excited state properties by solving for the linear-response function.

\[\delta n(r,\omega) = \int \chi(r,r',\omega)\,\delta v_{appl}(r',\omega)\,\,dr'\]

Lehman reresentation of the response function, .. math:: chi(r,r’,omega) = sum_{n}frac{n_{0n}^*(r)n_{0n}(r’)}{omega-Omega_{n}+ieta}-frac{n_{0n}^*(r)n_{0n}(r’)}{omega+Omega_{n}+ieta}

where \(n_{0n}(r) = \langle \Psi_0 | \hat n(r) | \Psi_n \rangle\) is the ground-to-$n^text{th}$ excited state transition density, and \(\Omega_n = E_n - E_0\).

Applying the Runge-Gross theorem,

\[\delta n(r,\omega) = \int \chi(r,r',\omega)\,\delta v_{appl}(r',\omega)\,\,dr' = \int \chi_s(r,r',\omega)\,\delta v_{s}(r',\omega)\,\,dr' = \int \chi_B(r,r',\omega)\,\delta v_{B}(r',\omega)\,\,dr'\]

where the variation of the KS potential can be decomposed into:

\[\delta v_{s}(r,\omega) = \delta v_{H}[n](r,\omega) + \delta v_{xc}[n](r,\omega) + \delta v_{appl}(r,\omega)\]

and the variation of the bosonic potential can be decomposed into crucial contributions

\[\delta v_{B}(r,\omega) = \delta v_s[n](r,\omega) + \delta v_P[n](r,\omega)\]

which can be derived from their dependence on the density

\[\delta v_{H}(r,\omega) = \int \frac{\delta v_{H}(r,\omega)}{\delta n(r',\omega)} \, \delta n(r',\omega) \, dr' = \int \frac{1}{|r-r'|}\, \delta n(r',\omega)\,dr'\]
\[\delta v_{xc}(r,\omega) = \int \frac{\delta v_{xc}(r,\omega)}{\delta n(r',\omega)} \, \delta n(r',\omega) \, dr' = \int f_{xc}(r,r',\omega)\, \delta n(r',\omega)\,dr'\]
\[\delta v_{P}(r,\omega) = \int \frac{\delta v_{P}(r,\omega)}{\delta n(r',\omega)} \, \delta n(r',\omega) \, dr' = \int f_{P}(r,r',\omega)\, \delta n(r',\omega)\,dr'\]

The Dyson equation for the response function between the real system and the KS system then given by:

\[-\chi^{-1} = -\chi_s^{-1} + \frac{1}{|r-r'|} + f_{xc}(r,r',\omega)\]

while the Dyson equation for the response function between the real system and the OF system then given by:

\[-\chi^{-1} = -\chi_B^{-1} + \frac{1}{|r-r'|} + f_{xc}(r,r',\omega) + f_P(r,r',\omega)\]

The Dyson equation turns into a matrix equation, known as the Casida equation, given by:

\[\begin{split}\left( \begin{array}{cc} A_{iajb}(\omega) & B_{iajb}(\omega) \\ B_{iajb}(\omega) & A_{iajb}(\omega) \end{array} \right) \left( \begin{array}{c} X_{ia} \\ Y_{ia} \end{array} \right) = \omega \left( \begin{array}{cc} 1 & 0 \\ 0 & -1 \end{array} \right) \left( \begin{array}{c} X_{ia} \\ Y_{ia} \end{array} \right)\end{split}\]

where defining \(\omega_{ia} = \varepsilon_a - \varepsilon_i\),

\[A_{iajb}(\omega) = \delta_{ij}\delta_{ab} \omega_{ia} + \big\langle \phi_i \phi_a \big| K(r,r',\omega) \big| \phi_j \phi_b \big\rangle\]

and

\[B_{iajb}(\omega) = \big\langle \phi_i \phi_a \big| K(r,r',\omega) \big| \phi_j \phi_b \big\rangle\]

where in KS \(K(r,r',\omega)\) is the kernel given by:

\[K(r,r',\omega) = \frac{1}{|r-r'|} + f_{xc}(r,r',\omega)\]

while in OF \(K(r,r',\omega)\) is the kernel given by:

\[K(r,r',\omega) = \frac{1}{|r-r'|} + f_{xc}(r,r',\omega) + f_P(r,r',\omega)\]

Collective modes such as plasmons are built from many particle–hole pairs in the fermionic KS picture but appear more directly in the bosonic reference, which can make Casida matrices much smaller for plasmon-dominated spectra while describing the same interacting response once approximations are controlled. Approximate \(f_P\) (e.g. Thomas–Fermi–von Weizsäcker in linear response) shifts peaks slightly compared with KS-TDDFT; Landau damping needs nonadiabatic kernels—DFTpy offers the nonadiabatic Pauli (JP) formulation for that regime.

In frequency space, DFTpy performs Casida TD-OFDFT using the modules Casida, CasidaRunner, and Hamiltonian on top of a converged ground-state density (see the linear-response tutorial lr-ofdft-tutorial.ipynb); this parallels real-time propagation in td-ofdft-tutorial.ipynb.

Short note on the implementation

In DFTpy, the electron density is represented in a discrete set of points given by a Cartesian grid and contained in a simulation cell that is specified by 3 lattice vectors. The number of grid points and the cell size are regulated by the user. The Cartesian grid allows for an efficient parallelization of data and work (we use mpi4py), and for the exploitation of Fast Fourier Transforms for solving convolution integrals (such as the one needed to compute \(E_H[n]\)). Either NumPy.fft or PyFFT are used depending on user input.

References