Jump to content

Requests for technical support from the VASP team should be posted in the VASP Forum.

Category:Time-dependent density functional theory: Difference between revisions

From VASP Wiki
Tal (talk | contribs)
Move Scaling into Casida TDDFT as L3 subsection; add explicit operation counts for building and solving
Tal (talk | contribs)
Replace Exchange-correlation kernel section with Time-evolution/real-time TDDFT section; remove scaling table
Line 60: Line 60:
* '''Exact diagonalization''' ({{TAG|IBSE}}{{=}}2): scales as <math>N_{\rm rank}^3</math>, or <math>N^6</math> with the system size.
* '''Exact diagonalization''' ({{TAG|IBSE}}{{=}}2): scales as <math>N_{\rm rank}^3</math>, or <math>N^6</math> with the system size.
* '''Time evolution''' ({{TAG|IBSE}}{{=}}1) and '''Lanczos''' ({{TAG|IBSE}}{{=}}3): scale as <math>N_{\rm rank}^2 \times m</math>, or <math>N^4</math> with the system size, where <math>m</math> is the number of time steps or iterations. The value of <math>m</math> depends on the required precision ({{TAG|BSEPREC}}) and the broadening ({{TAG|CSHIFT}}).
* '''Time evolution''' ({{TAG|IBSE}}{{=}}1) and '''Lanczos''' ({{TAG|IBSE}}{{=}}3): scale as <math>N_{\rm rank}^2 \times m</math>, or <math>N^4</math> with the system size, where <math>m</math> is the number of time steps or iterations. The value of <math>m</math> depends on the required precision ({{TAG|BSEPREC}}) and the broadening ({{TAG|CSHIFT}}).
{| class="wikitable"
|-
! !! Time evolution / Lanczos !! Exact diagonalization
|-
| Memory || <math>N_{\mathbf{k}} \times (N_v + N_c) \times N_G</math> || <math>(N_{\mathbf{k}} \times N_v \times N_c)^2</math>
|-
| Compute time || <math>N_{\mathbf{k}} \times N_v \times N_G</math> || <math>(N_{\mathbf{k}} \times N_v \times N_c)^3</math>
|-
| || <math>+ N_{\mathbf{k}} \times N_v \times N_c \times N_G</math> ||
|-
| Nonlocal exchange || <math>+ N_{\mathbf{k}}^2 \times N_v^2 \times N_G</math> || ⋯
|}


The following features are currently supported:
The following features are currently supported:
Line 79: Line 66:
* Calculations with [[hybrid functional|hybrid functionals]] and range-separated hybrids
* Calculations with [[hybrid functional|hybrid functionals]] and range-separated hybrids


== Exchange-correlation kernel ==
== Time-evolution or real-time TDDFT ==
 
The exchange-correlation kernel <math>f_\mathrm{xc}</math> determines how electron-hole interactions beyond the [[random phase approximation]] are described in TDDFT. The choice of kernel is tied to the exchange-correlation functional used in the ground-state calculation and is controlled by the tags {{TAG|LHARTREE}}, {{TAG|LADDER}}, and {{TAG|LFXC}}.
 
* '''Local and semilocal kernels (ALDA, APBE)''': obtained as the second functional derivative of an LDA or PBE exchange-correlation energy. These kernels are computationally cheap and work well for plasmons and metallic systems, but lack the long-range <math>-1/q^2</math> behavior and therefore fail to describe bound excitons in semiconductors and insulators.
 
* '''Hybrid-functional kernels''': when a fraction of exact exchange is included in the ground-state functional (e.g., PBE0 or HSE), the corresponding TDDFT kernel inherits a long-range non-local exchange contribution. This restores the <math>-1/q^2</math> behavior and allows for an approximate description of excitonic effects. The fraction of exact exchange is controlled by {{TAG|AEXX}}, and the range-separation parameter by {{TAG|HFSCREEN}}. When a hybrid functional is used, {{TAG|LADDER}} must be set to .TRUE. so that the non-local exchange contribution of the kernel (the ladder diagrams) is actually included in the time propagation; otherwise the calculation only contains the local part of <math>f_\mathrm{xc}</math> and the excitonic effects from the hybrid are lost.
 
* '''Screened exchange (ladder diagrams)''': enabling {{TAG|LADDER}} adds the screened Coulomb interaction <math>W</math>, yielding the proper long-range electron-hole attraction. In the time-evolution implementation, <math>W</math> is currently obtained from a [[Improving the dielectric function#Model-BSE|model dielectric function]] controlled by {{TAG|LMODELHF}}, {{TAG|AEXX}}, and {{TAG|HFSCREEN}}.


The combination of these three switches selects the approximation level: setting all three to .FALSE. yields the [[independent-particle approximation]] (equivalent to {{TAG|LOPTICS}}); enabling only {{TAG|LHARTREE}} gives the [[random phase approximation|RPA]]; further enabling {{TAG|LFXC}} or {{TAG|LADDER}} adds the beyond-RPA contributions.
An alternative to solving the Casida equation is to compute the frequency-dependent response via real-time propagation of the Kohn-Sham orbitals ({{TAG|ALGO}}{{=}}TIMEEV){{cite|sander:jcp:2017}}. The starting point is the time-dependent Kohn-Sham equation,
::<math>
\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>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>
\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>\mu_{cv\mathbf{k}}^j</math> is the dipole matrix element and <math>c^i_{cv\mathbf k}(t)</math> are the time-dependent expansion coefficients. This approach avoids building and storing the full excitonic Hamiltonian and scales as <math>N_{\rm rank}^2</math>, making it the method of choice for large systems. The detailed theory and propagation algorithm are described on the [[Construction:Time-dependent density-functional theory|theory page]].


== Additional resources ==
== Additional resources ==

Revision as of 12:11, 15 June 2026

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 implements TDDFT primarily through the Casida-equation formulation (ALGO=TDHF), which builds an excitonic Hamiltonian and solves it using exact diagonalization, time-evolution, or Lanczos algorithms. The detailed theoretical background is given in the theory page.

Casida TDDFT

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 with excitonic effects. 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 excitonic 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)

The Casida equation is solved via real-time propagation. A Dirac delta pulse of the electric field simultaneously excites all valence-to-conduction transitions, and the time-dependent dipole moments are propagated forward in time. 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 excitonic 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:

Time-evolution or real-time TDDFT

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 excitonic Hamiltonian and scales as [math]\displaystyle{ N_{\rm rank}^2 }[/math], making it the method of choice for large systems. The detailed theory and propagation algorithm are described on the theory page.

Additional resources

Lectures

Tutorials

How to

References

Pages in category "Time-dependent density functional theory"

The following 3 pages are in this category, out of 3 total.