How does DFTpy work?

Orbital-Free DFT

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]\). Namely,

\[n_0 = \text{argmin}_n \left[ E[n] - \mu \left( \int n(\mathbf{r}) d\mathbf{r} - N \right) \right]\]

where \(N\) is the number of valence electrons in the system. Both \(n_0(\mathbf{r})\) and \(\mu\) are determined during the minimization. The ground state energy is, \(E_0 = E[n_0]\).

In practice, the above minimization can only be carried out if the ground state energy functional is known as a pure functional of the density. The energy functional is a sum of several terms:

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

where

  • \(T_s[n]\): noninteracting kinetic energy or KEDF.

  • \(E_{xc}[n]\): exchange-correlation energy or EXC.

  • \(E_{H}[n]=\frac{1}{2}\int \frac{n(\mathbf{r})n(\mathbf{r}^\prime)}{|\mathbf{r}-\mathbf{r}^\prime|}d\mathbf{r} d\mathbf{r}^\prime\): Hartree energy.

  • \(v_{ext}(\mathbf{r})\): the external potential (typically the electron-ion interaction).

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 Orbital-Free DFT

DFTpy can also describe systems out of equilibrium by propagating them in real time as well as in frequency space finding the roots of the frequency dependent polarizability (Casida). Because of the OF-DFT foundation, DFTpy implements the so-called time-dependent orbital-free DFT (td-OF-DFT) approach whereby a single Bosonic wavefunction, \(\Psi(\mathbf{r},t)\) is propagated with a time-dependent KS-like Hamiltonian. Namely,

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

where

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

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

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

where \(v_s(\mathbf{r},t)=v_{ext}(\mathbf{r},t)+v_H[n(t)](\mathbf{r},t)+v_{xc}[n(t)](\mathbf{r},t)\) where the adiabatic approximation has been invoked. The Pauli potential is given by adiabatic and nonadiabatic contributions,

\[v_P(\mathbf{r},t)=v_P^{ad}(\mathbf{r},t)+v_P^{nad}(\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 litarature. 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.

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