Category:Time-dependent density functional theory

Time-dependent density-functional theory (TDDFT) extends density-functional theory to time-varying external potentials, enabling the computation of neutral electronic excitations and frequency-dependent response functions. VASP provides multiple implementations of TDDFT, each with its own advantages and disadvantages, so that the best algorithm can be selected based on the problem at hand. The detailed theoretical background is given on the theory page.

Casida TDDFT (ALGO=TDHF)

The Casida formulation of TDDFT (ALGO=TDHF) recasts the linear-response problem as a non-Hermitian eigenvalue problem

[math]\displaystyle{ \left(\begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \end{array}\right)\left(\begin{array}{l} \mathbf{X}_\lambda \\ \mathbf{Y}_\lambda \end{array}\right)=\omega_\lambda\left(\begin{array}{cc} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & -\mathbf{1} \end{array}\right)\left(\begin{array}{l} \mathbf{X}_\lambda \\ \mathbf{Y}_\lambda \end{array}\right), }[/math]

with the same mathematical structure as the Bethe-Salpeter equation (BSE). The eigenvalues [math]\displaystyle{ \omega_\lambda }[/math] are the excitation energies, and the eigenvectors [math]\displaystyle{ \mathbf{X}_\lambda, \mathbf{Y}_\lambda }[/math] determine the oscillator strengths and the dielectric function. The difference with respect to BSE is that the beyond-RPA part of the kernel is described by the exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc} }[/math] instead of the screened Coulomb interaction [math]\displaystyle{ W }[/math].

Casida TDDFT can be performed using DFT or hybrid-functional orbitals and eigenvalues. VASP provides three algorithms for solving the Casida equation, selected via IBSE:

Exact diagonalization (IBSE=2)

The Hamiltonian is diagonalized exactly. The excitation energies and oscillator strengths are obtained directly from the eigenvalues and eigenvectors, which makes this approach particularly useful for analyzing individual excitons. The macroscopic dielectric function is obtained from the spectral representation:

[math]\displaystyle{ \varepsilon_M(\omega) = 1 + \frac{4\pi}{\Omega} \sum_\lambda \left|\sum_{cv\mathbf k} \mu_{cv\mathbf k} X_\lambda^{cv\mathbf k}\right|^2 \left[\frac{1}{\omega + \omega_\lambda + \mathrm i\eta} - \frac{1}{\omega - \omega_\lambda + \mathrm i\eta}\right], }[/math]

where [math]\displaystyle{ \mu_{cv\mathbf{k}}^j=\frac{\langle c\mathbf{k}|v_j|v\mathbf{k}\rangle}{\varepsilon_c(\mathbf{k})-\varepsilon_v(\mathbf{k})} }[/math] is the dipole matrix element, [math]\displaystyle{ X_\lambda^{cv\mathbf k} }[/math] are the eigenvector components, and [math]\displaystyle{ \omega_\lambda }[/math] are the excitation energies.

Time-evolution algorithm (IBSE=1)

This algorithm is based on the same time-evolution approach described in the real-time TDDFT section, but applied within the Casida framework. The key difference is that the Hamiltonian in transition space is built once and never updated during the propagation. Only the time-dependent dipole vector [math]\displaystyle{ |\mu_{cv\mathbf k}(t)\rangle }[/math] is propagated forward in time using the fixed Hamiltonian.

A Dirac delta pulse of the electric field simultaneously excites all valence-to-conduction transitions, and the dielectric function is found via a Fourier transform^[1]:

[math]\displaystyle{ \varepsilon_M(\omega)=1-\frac{4\pi}{\Omega}\int_0^{\infty} \mathrm{d} t \sum_{c,v,\mathbf{k}}\left(\langle\mu_{cv\mathbf{k}}| \xi_{cv\mathbf{k}}(t)\rangle + \mathrm{c.c.}\right) e^{-\mathrm i(\omega-\mathrm i \eta) t}, }[/math]

where [math]\displaystyle{ \mu_{cv\mathbf k} }[/math] are the dipole moments and [math]\displaystyle{ |\xi_{cv\mathbf k}(t)\rangle }[/math] is the time-evolved dipole vector. The solution is strictly equivalent to that of the exact diagonalization for the dielectric function, but does not yield eigenvectors and so cannot be used directly for exciton analysis. Its main advantage is the quadratic scaling with [math]\displaystyle{ N_{\rm rank} }[/math]. The required number of propagation steps is controlled by the broadening CSHIFT and the maximum energy OMEGAMAX, and does not depend on the size of the Hamiltonian.

Lanczos algorithm (IBSE=3)

The dielectric function is expressed as a continued fraction

[math]\displaystyle{ \varepsilon_M(\omega) = 1 - \frac{4\pi}{\Omega}\cfrac{|u_0|^2}{(\omega - a_1 + \mathrm i\eta) - \cfrac{b_1^2}{(\omega -a_2 + \mathrm i\eta) - \cfrac{b_2^2}{...}}}, }[/math]

where [math]\displaystyle{ |u_0\rangle }[/math] is an initial guess vector computed from the dipole moments. The [math]\displaystyle{ a }[/math] and [math]\displaystyle{ b }[/math] coefficients are evaluated iteratively, with the algorithm stopping once the difference between [math]\displaystyle{ \varepsilon(\omega) }[/math] from two consecutive iterations is below a threshold selected by BSEPREC. Because the starting vector is built from dipole moments, the Lanczos algorithm is sensitive only to optically active transitions and can reach convergence faster than other methods for larger matrices.

Scaling

As a first step, the algorithm builds the Hamiltonian of rank

[math]\displaystyle{ N_{\rm rank} = N_k \times N_c \times N_v, }[/math]

where [math]\displaystyle{ N_k }[/math], [math]\displaystyle{ N_c }[/math], and [math]\displaystyle{ N_v }[/math] are the number of k points, conduction bands, and valence bands. Building the Hamiltonian requires

[math]\displaystyle{ N_k \times N_q \times N_v^2 \times N_c^2 \times N_G }[/math]

operations, which scales as [math]\displaystyle{ N^4 }[/math]–[math]\displaystyle{ N^5 }[/math] with the system size.

In the second step, the equation is solved. The scaling depends on the algorithm:

• Exact diagonalization (IBSE=2): scales as [math]\displaystyle{ N_{\rm rank}^3 }[/math], or [math]\displaystyle{ N^6 }[/math] with the system size.
• Time evolution (IBSE=1) and Lanczos (IBSE=3): scale as [math]\displaystyle{ N_{\rm rank}^2 \times m }[/math], or [math]\displaystyle{ N^4 }[/math] with the system size, where [math]\displaystyle{ m }[/math] is the number of time steps or iterations. The value of [math]\displaystyle{ m }[/math] depends on the required precision (BSEPREC) and the broadening (CSHIFT).

The following features are currently supported:

• Calculating the dielectric function and eigenvectors (IBSE=2)
• Tamm-Dancoff approximation
• Calculations with hybrid functionals and range-separated hybrids

Time-evolution or real-time TDDFT (ALGO=TIMEEV)

An alternative to solving the Casida equation is to compute the frequency-dependent response via real-time propagation of the Kohn-Sham orbitals (ALGO=TIMEEV)^[1]. The starting point is the time-dependent Kohn-Sham equation,

[math]\displaystyle{ \mathrm i \frac{\partial}{\partial t}\left|\phi_{v\mathbf k}[n(t)]\right\rangle = \left[-\frac{\nabla^2}{2} + V_{\mathrm H}[n(t)] + V_{\mathrm{xc}}[n(t)] + V_\mathrm{ext}(t)\right]\left|\phi_{v\mathbf k}[n(t)]\right\rangle, }[/math]

where all potentials are functionals of the time-evolving density [math]\displaystyle{ n(\mathbf r, t) }[/math]. A delta-like perturbation is applied to the ground-state system, and the time-dependent coefficients are propagated forward in time. The dielectric function is then obtained from the Fourier transform of the time-dependent dipole moments:

[math]\displaystyle{ \varepsilon_{ij}(\omega) = \delta_{ij} - \frac{4\pi e^2}{\Omega} \int_0^\infty \mathrm d t \sum_{cv\mathbf k} \left(\langle\mu^j_{cv\mathbf k}|c^i_{cv\mathbf k}(t)\rangle + \mathrm{c.c.}\right) e^{-\mathrm i(\omega - \mathrm i\eta)t}, }[/math]

where [math]\displaystyle{ \mu_{cv\mathbf{k}}^j }[/math] is the dipole matrix element and [math]\displaystyle{ c^i_{cv\mathbf k}(t) }[/math] are the time-dependent expansion coefficients. This approach avoids building and storing the full Hamiltonian, making it the method of choice for large systems. The detailed theory and propagation algorithm are described on the theory page.

Scaling

The memory requirement scales as [math]\displaystyle{ N_{\mathbf{k}} \times (N_v + N_c) \times N_G }[/math], since only the wavefunctions and their time-dependent coefficients need to be stored. The compute time per time step scales as

[math]\displaystyle{ N_{\mathbf{k}} \times N_v \times N_G + N_{\mathbf{k}} \times N_v \times N_c \times N_G, }[/math]

where the first term accounts for updating the Hamiltonian and the second for propagating the coefficients. When nonlocal exchange is included (LADDER=.TRUE.), an additional [math]\displaystyle{ N_{\mathbf{k}}^2 \times N_v^2 \times N_G }[/math] contribution arises from evaluating the exchange integrals at each time step.

Dyson-equation TDDFT (ALGO=CHI)

Instead of recasting the problem as an eigenvalue equation, TDDFT can also be solved by directly evaluating the two-point Dyson equation for the density-density response function [math]\displaystyle{ \chi }[/math] (ALGO=CHI). The macroscopic dielectric function is then obtained from

[math]\displaystyle{ \varepsilon_M(\omega) = 1 - v(\mathbf{G}=0) \, \chi_{\mathbf{G}=0,\mathbf{G}'=0}(\mathbf{q},\omega), }[/math]

where [math]\displaystyle{ v(\mathbf{G}) }[/math] is the Coulomb kernel and [math]\displaystyle{ \chi_{\mathbf{G},\mathbf{G}'}(\mathbf{q},\omega) }[/math] is the interacting response function in reciprocal space. This approach gives access to the full matrix [math]\displaystyle{ \chi_{\mathbf{G},\mathbf{G}'}(\mathbf{q},\omega) }[/math] and is used, for instance, to compute the screened Coulomb interaction [math]\displaystyle{ W }[/math] needed in GW calculations. However, the Dyson equation must be inverted at every frequency point, which makes it expensive for computing spectra with fine resolution. Since this approach does not invoke the Tamm-Dancoff approximation, the RPA dielectric function obtained this way is equivalent to the Casida RPA dielectric function. The only nonlocal exchange-correlation kernel currently supported is the nanoquanta kernel (LFXC).

Scaling

The compute time scales as [math]\displaystyle{ N^4 }[/math] with the system size. At each frequency point, the independent-particle response [math]\displaystyle{ \chi_0(\mathbf{q},\omega) }[/math] is constructed from the sum over transitions, which scales as [math]\displaystyle{ N_{\mathbf{k}} \times N_v \times N_c \times N_G^2 }[/math]. The subsequent matrix inversion of the Dyson equation scales as [math]\displaystyle{ N_G^3 }[/math]. Because these operations must be repeated for every frequency point, the total cost grows linearly with the number of frequency points.

Additional resources

Lectures

• Lecture on TDDFT theory and calculations.

Tutorials

• Tutorial calculating optical absorption of C diamond using TDDDH.
• Tutorial on the efficient Brillouin zone sampling using TDDDH and exciton analysis using TDDDH.

How to

• Practical guide for solving the Casida equation via diagonalization: TDDFT calculations.
• Practical guide for real-time TDDFT calculations: Time-evolution algorithm.

References