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)
No edit summary
Tal (talk | contribs)
 
(52 intermediate revisions by 3 users not shown)
Line 1: Line 1:
The formalism of the Bethe-Salpeter equation (BSE) allows for calculating the polarizability with the electron-hole interaction and constitutes the state of the art for calculating absorption spectra in solids.
'''Time-dependent density-functional theory''' (TDDFT) extends [[:Category:Exchange-correlation functionals|density-functional theory]] to time-varying external potentials, enabling the computation of neutral electronic excitations and frequency-dependent response functions. The accuracy of TDDFT strongly depends on the choice of the exchange-correlation kernel <math>f_\mathrm{xc}</math> ({{TAG|LFXC}}, {{TAG|LADDER}}) and if the nonlocal limit of <math>f_\mathrm{xc}</math> is properly treated, excitonic effects can be captured with good accuracy{{cite|tal:prr:2020}}. 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 [[Time-dependent density-functional theory|theory page]].


== Theory ==
* Lecture on {{Video|bse:alexey:2026|TDDFT theory and calculations}}.


Time-dependent density-functional theory (TDDFT) is an extension of DFT to address excited-state properties, dynamics, and spectroscopy. In principle, TDDFT is an exact theory for neutral electronic excitations, however, similarly to DFT, the exchange-correlation functional is unknown and needs to be approximated.
== Casida TDDFT ({{TAG|ALGO|TDHF}}) ==
 
In the linear response approximation, we split the external potential into a static term and a time-dependent perturbation <math display="block">v(r,t) =v(r)+\delta v(r,t),</math> where the perturbation term is much smaller than the static potential <math display="inline">\delta
v(r,t) \ll v(r)</math>. In this case the Hohenberg–Kohn and Runge–Gross theorems state the correspondence <math display="inline">\delta \rho(r,t) \Leftrightarrow \delta v(r,t)</math>. A TDDFT calculation is a two-step procedure: first, we perform an ordinary DFT calculation with a static external potential <math display="inline">v(r)</math> and then we perform a TDDFT calculation of the density variation <math display="inline">\delta \rho(r,t)</math> corresponding to the external time-dependent perturbation <math display="inline">\delta v(r,t)</math>. From <math display="inline">\delta \rho(r,t)</math> we can calculate the polarizability of the system <math display="inline">\chi</math> using <math display="block">\delta \rho(r_1,t_1)= \int dr_2dt_2 \chi(r_1,t_1,r_2,t_2)\delta v(r_2,t_2).</math>
 
Following a Kohn-Sham (KS) scheme we assume that the density response of KS system is equivalent to that of the real system, i.e., <math display="inline">\delta \rho = \delta \rho^{\rm
KS}</math>, in response to an effective KS petrubation
 
<math display="block">\delta v^{\mathrm{KS}}(x)=\delta v(x)+\delta v_{\mathrm{H}}(x)+\delta
v_{\mathrm{xc}}(x).</math> Here, <math display="inline">\delta v(r_1,t_1)</math> is the real external perturbation, the Hartree term <math display="inline">\delta v_H(x)</math> is <math display="block">\delta v_{\mathrm{H}}(x)=\int dr_2dt_2 V(r_1,t_1,r_2,t_2)\delta \rho(r_2,t_2)</math> and the exchange-correlation term <math display="block">\delta v_{\mathrm{xc}}(x)=\int dr_2dt_2 f_{\rm xc}[\rho](r_1,t_1,r_2,t_2)\delta \rho(r_2,t_2).</math> The main challenge lays in finding an accurate approximation for the exchange-correlation kernel <math display="block">f_{\rm xc}[\rho](r_1,t_1,r_2,t_2)=\frac{\delta v_{\rm xc}[\rho](r_1,t_1)}{\delta \rho(r_2,t_2)}</math>
 
The response of the non-interacting KS particles is then <math display="block">\delta \rho(r_1,t_1)= \int dr_2dt_2 \chi^{\rm KS}(r_1,t_1,r_2,t_2)\delta v^{\rm KS}(r_2,t_2).</math>
 
Then, writing the Adler-Wiser expression in reciprocal space and frequency domain we can find the response function of the KS system
 
<math display="block">\begin{gathered}
\chi_{\mathbf{G}, \mathbf{G}^{\prime}}^{K S}(\mathbf{q}, \omega)
=-\frac{1}{V} \sum_{n \mathbf{k}} \sum_{m \mathbf{k}^{\prime}} 2 f_{n
\mathbf{k}}\left(1-f_{m \mathbf{k}^{\prime}}\right)\left(\frac{\left\langle m
\mathbf{k}^{\prime}\left|e^{i(\mathbf{q}+\mathbf{G}) \mathbf{r}}\right| n
\mathbf{k}\right\rangle\left\langle n
\mathbf{k}\left|e^{-i\left(\mathbf{q}+\mathbf{G}^{\prime}\right)
\mathbf{r}^{\prime}}\right| m \mathbf{k}^{\prime}\right\rangle}{\epsilon_{m
\mathbf{k}^{\prime}}-\epsilon_{n \mathbf{k}}-\bar{\omega}}+\right. \\
\left.+\frac{\left\langle n \mathbf{k}\left|e^{i(\mathbf{q}+\mathbf{G})
\mathbf{r}}\right| m \mathbf{k}^{\prime}\right\rangle\left\langle m
\mathbf{k}^{\prime}\left|e^{-i\left(\mathbf{q}+\mathbf{G}^{\prime}\right)
\mathbf{r}^{\prime}}\right| n \mathbf{k}\right\rangle}{\epsilon_{m
\mathbf{k}^{\prime}}-\epsilon_{n \mathbf{k}}+\bar{\omega}}\right) \\
\end{gathered}</math>
 
and using the Dyson equation for the polarizability we find the polarizability of the real system
 
<math display="block">\chi = \chi^{\rm KS}+ \chi^{\rm KS}(v+f_{\rm xc})\chi.</math>
 
The exciation frequencies of the system can be extracted from the analytic structure of polarizability <math display="inline">\chi</math>.
 
Finally, the dielectric function is found <math display="block">\varepsilon^{-1}_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega)=
\delta_{\mathbf{G}, \mathbf{G}^{\prime}}+\frac{4 \pi e^2}{|\mathbf{G}+\mathbf{q}|\left|\mathbf{G}^{\prime}+\mathbf{q}\right|} \chi_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega)</math>
 
# Casida equation Alternatively, the excitation energies <math display="inline">\omega_\lambda</math> of the real system can be found by mapping Eq. (1) onto an eigenvalue problem
 
<math display="block">\left(\begin{array}{cc}
A & B \\
B^* & A^*
\end{array}\right)\left(\begin{array}{l}
X_\lambda \\
Y_\lambda
\end{array}\right)=\omega_\lambda\left(\begin{array}{cc}
1 & 0 \\
0 & -1
\end{array}\right)\left(\begin{array}{l}
X_\lambda \\
Y_\lambda
\end{array}\right)</math>
 
The structure of the Casida equation is very similar to that of the Bethe-Salpeter equation. And similarly to BSE, the standard way to solve the Casida equation is to to neglect the coupling terms <math display="inline">B</math> and <math display="inline">B^*</math>, i.e., the Tamm-Dancoff approximation
 
<math display="block">AX_\lambda=\omega_\lambda X_\lambda~.</math>
 
The interaction between in TDDFT is described by the bare Coulomb <math display="inline">V_\mathbf{G}</math> and the exchange-correlation kernel
 
<math display="block">\begin{gathered}
A_{vc\mathbf{k}}^{v'c'\mathbf{k}'} =
(\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'}\delta_{\mathbf{kk}'}+
\frac{2}{\Omega}\sum_{\mathbf{G}\neq0}\bar{V}_\mathbf{G}(\mathbf{q})\langle
c\mathbf{k}|e^{i\mathbf{Gr}}|v\mathbf{k}\rangle\langle
v'\mathbf{k}'|e^{-i\mathbf{Gr}}|c'\mathbf{k}'\rangle\\
+
\frac{2}{\Omega}\sum_{\mathbf{G}\neq0}f^{\rm xc}_{\mathbf{G,G'}}\langle
c\mathbf{k}|e^{i\mathbf{Gr}}|v\mathbf{k}\rangle\langle
v'\mathbf{k}'|e^{-i\mathbf{Gr}}|c'\mathbf{k}'\rangle.
\end{gathered}</math>
 
If xc potential includes the non-local exact exchange contribution, an additional term will appear
 
<math display="block">\begin{gathered}
A_{vc\mathbf{k}}^{v'c'\mathbf{k}'} =
(\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'}\delta_{\mathbf{kk}'}+
\frac{2}{\Omega}\sum_{\mathbf{G}\neq0}\bar{V}_\mathbf{G}(\mathbf{q})\langle
c\mathbf{k}|e^{i\mathbf{Gr}}|v\mathbf{k}\rangle\langle
v'\mathbf{k}'|e^{-i\mathbf{Gr}}|c'\mathbf{k}'\rangle\\
+ \frac{2}{\Omega}\sum_{\mathbf{G}\neq0}f^{\rm xc,loc}_{\mathbf{G,G'}}\langle
c\mathbf{k}|e^{i\mathbf{Gr}}|v\mathbf{k}\rangle\langle
v'\mathbf{k}'|e^{-i\mathbf{Gr}}|c'\mathbf{k}'\rangle
\\-\frac{2}{\Omega}\sum_{\mathbf{G,G}'}c_{\rm x}(\mathbf{q+G})V_{\mathbf{G}}(\mathbf{q})\delta_{\mathbf{q,k-k}'}
\langle c\mathbf{k}|e^{i(\mathbf{q+G})}|c'\mathbf{k}'\rangle \langle
v'\mathbf{k}'|e^{-i(\mathbf{q+G})}|v\mathbf{k}\rangle.
\end{gathered}</math>
 
Here, <math display="inline">c_x(\mathbf{q+G})</math> is the range-dependent fraction of the exact exchange potential.
 
# Common approximations Neglecting both interaction terms, i.e., <math display="inline">v</math> and <math display="inline">f_{\rm xc}</math> yields the independent particle approximation.
 
=== Bethe-Salpeter equation ===
 
In the BSE, the excitation energies correspond to the eigenvalues <math>\omega_\lambda</math> of the following linear problem


The Casida formulation of TDDFT recasts the linear-response problem as a non-Hermitian eigenvalue problem
::<math>
::<math>
\left(\begin{array}{cc}
\left(\begin{array}{cc}
Line 114: Line 19:
\mathbf{X}_\lambda \\
\mathbf{X}_\lambda \\
\mathbf{Y}_\lambda
\mathbf{Y}_\lambda
\end{array}\right)~.
\end{array}\right),
</math>
</math>
with the same mathematical structure as the [[Bethe-Salpeter equation]] (BSE). The eigenvalues <math>\omega_\lambda</math> are the excitation energies, and the eigenvectors <math>\mathbf{X}_\lambda, \mathbf{Y}_\lambda</math> determine the oscillator strengths and the dielectric function.
VASP provides three algorithms for solving the Casida equation, selected via {{TAG|IBSE}}:


=== Exact diagonalization ({{TAG|IBSE|2}}) ===


The matrices <math>A</math> and <math>A^*</math> describe the resonant and anti-resonant transitions between the occupied <math>v,v'</math> and unoccupied <math>c,c'</math> states
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>
::<math>
A_{vc}^{v'c'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'} + \langle cv'|V|vc'\rangle - \langle cv'|W|c'v\rangle.
\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>
</math>
The energies and orbitals of these states are usually obtained in a <math>G_0W_0</math> calculation, but DFT and Hybrid functional calculations can be used as well.
where <math>\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>X_\lambda^{cv\mathbf k}</math> are the eigenvector components, and <math>\omega_\lambda</math> are the excitation energies.
The electron-electron interaction and electron-hole interaction are described via the bare Coulomb <math>V</math> and the screened potential <math>W</math>.


The coupling between resonant and anti-resonant terms is described via terms <math>B</math> and <math>B^*</math>
=== Time-evolution algorithm ({{TAG|IBSE|1}}) ===
::<math>
B_{vc}^{v'c'} = \langle vv'|V|cc'\rangle - \langle vv'|W|c'c\rangle.
</math>
Due to the presence of this coupling, the Bethe-Salpeter Hamiltonian is non-Hermitian.


=== Tamm-Dancoff approximation ===
This algorithm is based on the same time-evolution approach described in the [[#Time-evolution or real-time TDDFT (ALGO=TIMEEV)|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>|\mu_{cv\mathbf k}(t)\rangle</math> is propagated forward in time using the fixed Hamiltonian.
A common approximation to the BSE is the Tamm-Dancoff approximation (TDA), which neglects the coupling between resonant and anti-resonant terms, i.e., <math>B</math> and <math>B^*</math>.  
Hence, the TDA reduces the BSE to a Hermitian problem


The dielectric function is found via a Fourier transform{{cite|sander:jcp:2017}}:
::<math>
::<math>
AX_\lambda=\omega_\lambda X_\lambda~.
\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>
</math>
where <math>\mu_{cv\mathbf k}</math> are the dipole moments and <math>|\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.


In reciprocal space, the matrix <math>A</math> is written as 
=== Lanczos algorithm ({{TAG|IBSE|3}}) ===


The dielectric function is expressed as a continued fraction
::<math>
::<math>
A_{vc\mathbf{k}}^{v'c'\mathbf{k}'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'}\delta_{\mathbf{kk}'}+
\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)  
\frac{2}{\Omega}\sum_{\mathbf{G}\neq0}\bar{V}_\mathbf{G}(\mathbf{q})\langle c\mathbf{k}|e^{i\mathbf{Gr}}|v\mathbf{k}\rangle\langle v'\mathbf{k}'|e^{-i\mathbf{Gr}}|c'\mathbf{k}'\rangle
    - \cfrac{b_2^2}{...}}},
  -\frac{2}{\Omega}\sum_{\mathbf{G,G}'}W_{\mathbf{G,G}'}(\mathbf{q},\omega)\delta_{\mathbf{q,k-k}'}
\langle c\mathbf{k}|e^{i(\mathbf{q+G})}|c'\mathbf{k}'\rangle
\langle v'\mathbf{k}'|e^{-i(\mathbf{q+G})}|v\mathbf{k}\rangle,
</math>
</math>
where <math>|u_0\rangle</math> is an initial guess vector computed from the dipole moments. The <math>a</math> and <math>b</math> coefficients are evaluated iteratively, with the algorithm stopping once the difference between <math>\varepsilon(\omega)</math> from two consecutive iterations is below a threshold selected by {{TAG|BSEPREC}}.


where <math>\Omega</math> is the cell volume, <math>\bar{V}</math> is the bare Coulomb potential without the long-range part
=== Scaling with the system size===


::<math>
Building the Hamiltonian scales as <math>N^4</math>–<math>N^5</math> with the system size. Solving the equation scales as <math>N^6</math> for exact diagonalization ({{TAG|IBSE|2}}) and <math>N^4</math> for the time-evolution ({{TAG|IBSE|1}}) and Lanczos ({{TAG|IBSE|3}}) algorithms.
\bar{V}_{\mathbf{G}}(\mathbf{q})=\begin{cases}
    0 & \text { if } G=0 \\
    V_{\mathbf{G}}(\mathbf{q})=\frac{4 \pi}{|\mathbf{q}+\mathbf{G}|^2} & \text { else }
\end{cases}~,
</math>
and the screened Coulomb potential
<math>
W_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega)=\frac{4 \pi \epsilon_{\mathbf{G}, \mathbf{G}^{\prime}}^{-1}(\mathbf{q}, \omega)}{|\mathbf{q}+\mathbf{G}|\left|\mathbf{q}+\mathbf{G}^{\prime}\right|}.
</math>


Here, the dielectric function <math>\epsilon_\mathbf{G,G'}(\mathbf{q})</math> describes the screening in <math>W</math> within the random-phase approximation (RPA)
== Time-evolution or real-time TDDFT ({{TAG|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 {{cite|sander:jcp:2017}}. The starting point is the time-dependent Kohn-Sham equation,
::<math>
::<math>
\epsilon_{\mathbf{G}, \mathbf{G}^{\prime}}^{-1}(\mathbf{q}, \omega)=\delta_{\mathbf{G}, \mathbf{G}^{\prime}}+\frac{4 \pi}{|\mathbf{q}+\mathbf{G}|^2} \chi_{\mathbf{G}, \mathbf{G}^{\prime}}^{\mathrm{RPA}}(\mathbf{q}, \omega).
\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>
</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:
Although the dielectric function is frequency-dependent, the static approximation <math>W_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega=0)</math> is considered a standard for practical BSE calculations.  
=== Macroscopic dielectric function ===
The macroscopic dielectric which accounts for the excitonic effects is found via eigenvalues <math>\omega_\lambda</math> and eigenvectors <math>X_\lambda</math> of the BSE
 
::<math>
::<math>
\epsilon_M(\mathbf{q}=0,\omega)=
\varepsilon_M(\omega) = 1 - \frac{4\pi e^2}{\Omega} \int_0^\infty \mathrm d t \sum_{cv\mathbf k} \left(\langle\mu_{cv\mathbf k}|c_{cv\mathbf k}(t)\rangle + \mathrm{c.c.}\right) e^{-\mathrm i(\omega - \mathrm i\eta)t},
1-\lim_{\mathbf{q}\rightarrow 0}v(q)\sum_{\lambda}
\left|\sum_{c,v,k}\langle c\mathbf{k}|e^{i\mathbf{qr}}|v\mathbf{k}\rangle X_\lambda^{cv\mathbf{k}}\right|^2
\times
\left(\frac{1}{\omega - \omega_\lambda + i\delta}\right)~.
</math>
</math>
where <math>\mu_{cv\mathbf{k}}</math> is the dipole matrix element and <math>c_{cv\mathbf k}(t)</math> are the time-dependent expansion coefficients. This approach avoids building and storing the full Hamiltonian in transition space. The detailed theory and propagation algorithm are described on the [[Time-dependent density-functional theory|theory page]].


== Scaling ==
=== Scaling with the system size===
The scaling of the BSE equation strongly limits its application for large systems. The main limiting factor is the diagonalization of the BSE Hamiltonian. The rank of the Hamiltonian is


:<math>N_{\rm rank} = N_k\times N_c\times N_v</math>,
The compute time per time step scales as <math>N^3</math> with the system size. When nonlocal exchange is included ({{TAG|LADDER}}{{=}}.TRUE.), an additional <math>N_{\mathbf{k}}^2 \times N_v^2 \times N_G</math> contribution arises from evaluating the exchange integrals at each time step.


where <math>N_k</math> is the number of k-points in the Brillouin zone and <math>N_c</math> and <math>N_v</math> are the number of conduction and valence bands, respectively.  The diagonalization of the matrix scales cubically with the matrix rank, i.e.,  <math>N_{\rm rank}^3</math>.
== Dyson-equation TDDFT ({{TAG|ALGO|CHI}}) ==


Despite the fact that this matrix diagonalization is usually the bottleneck for bigger systems, the construction of the BSE Hamiltonian also scales unfavorably and can play a dominant role in big systems, i.e.,
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>\chi</math>. The macroscopic dielectric function is then obtained from
::<math>
\varepsilon_M(\omega) = 1 - v(\mathbf{G}=0) \, \chi_{\mathbf{G}=0,\mathbf{G}'=0}(\mathbf{q},\omega),
</math>
where <math>v(\mathbf{G})</math> is the Coulomb kernel and <math>\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>\chi_{\mathbf{G},\mathbf{G}'}(\mathbf{q},\omega)</math> and is used, for instance, to compute the screened Coulomb interaction <math>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. The only nonlocal exchange-correlation kernel currently supported is the nanoquanta kernel ({{TAG|LFXC}}).


:<math>N_k\times N_q\times (N_v\times N_v\times N_G\times N_c\times N_c)</math>,
=== Scaling with the system size===


where <math>N_q</math> is the number of q-points and <math>N_G</math> number of G-vectors.
The compute time scales as <math>N^4</math> with the system size. The Dyson equation must be inverted at every frequency point ({{TAG|NOMEGA}}), so the total cost grows linearly with the number of frequency points.


== Solution of BSE ==
== Additional resources ==
=== Diagonalization ===
The exact diagonalization of the BSE Hamiltonian can be perform using various eigensolvers provided in ScaLAPACK, ELPA, and cuSolver libraries. The advantage of this approach is that the eigenvectors can be directly obtained and used for the analysis of the excitons.


The following features are currently supported:
=== Tutorials ===
* [[BSE calculations#Bethe-Salpeter equation calculation|Obtaining the spectra and eigenvectors]]
* Tutorial calculating optical absorption of C diamond using {{Tutorial|bse:e03|TDDDH}}.
* [[BSE calculations#Calculations beyond Tamm-Dancoff approximation|Calculations beyond Tamm-Dancoff approximation]]
* Tutorial on the {{Tutorial|bse:e09|efficient Brillouin zone sampling}} using TDDDH and {{Tutorial|bse:e10|exciton analysis}} using TDDDH.
* [[BSE calculations#Calculations at finite wavevectors|Calculations of <math>\varepsilon(\mathbf{q},\omega)</math> for <math>\mathbf{q}\neq0</math>]]
* Fatband plot


=== Time evolution ===
=== How to ===
The alternative approach is to formulate the BSE as the initial-value problem for the macroscopic polarizability. This approach converges to the same solution as the exact diagonalization and can be used for obtaining the absorption spectrum, but does not yield the eigenvectors, which can be limiting for the analysis of the excitons. The advantage of this approach is the quadratic scaling with the size of the BSE Hamiltonian <math>N_{rank}^2</math>.
* Practical guide for solving the Casida equation via diagonalization: [[TDDFT calculations]].
* Practical guide for real-time TDDFT calculations: [[Time-evolution algorithm]].


The following features are currently supported:
== References ==
* Obtaining the spectra
* Calculations beyond Tamm-Dancoff approximation
* Calculations of <math>\varepsilon(\mathbf{q},\omega)</math> for <math>\mathbf{q}\neq0</math>
 
== How to ==
 
* Practical guide for solving the BSE via diagonalization [[BSE calculations]]


== References ==
<references/>


[[Category:VASP|Bethe-Salpeter equation]][[Category:Many-body perturbation theory]]
[[Category:VASP|TDDFT]][[Category:Many-body perturbation theory]][[Category:Linear response]]

Latest revision as of 11:40, 17 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. The accuracy of TDDFT strongly depends on the choice of the exchange-correlation kernel [math]\displaystyle{ f_\mathrm{xc} }[/math] (LFXC, LADDER) and if the nonlocal limit of [math]\displaystyle{ f_\mathrm{xc} }[/math] is properly treated, excitonic effects can be captured with good accuracy[1]. 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 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.

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.

The dielectric function is found via a Fourier transform[2]:

[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.

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.

Scaling with the system size

Building the Hamiltonian scales as [math]\displaystyle{ N^4 }[/math][math]\displaystyle{ N^5 }[/math] with the system size. Solving the equation scales as [math]\displaystyle{ N^6 }[/math] for exact diagonalization (IBSE = 2) and [math]\displaystyle{ N^4 }[/math] for the time-evolution (IBSE = 1) and Lanczos (IBSE = 3) algorithms.

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 [2]. 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_M(\omega) = 1 - \frac{4\pi e^2}{\Omega} \int_0^\infty \mathrm d t \sum_{cv\mathbf k} \left(\langle\mu_{cv\mathbf k}|c_{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] is the dipole matrix element and [math]\displaystyle{ c_{cv\mathbf k}(t) }[/math] are the time-dependent expansion coefficients. This approach avoids building and storing the full Hamiltonian in transition space. The detailed theory and propagation algorithm are described on the theory page.

Scaling with the system size

The compute time per time step scales as [math]\displaystyle{ N^3 }[/math] with the system size. 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]. 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. The only nonlocal exchange-correlation kernel currently supported is the nanoquanta kernel (LFXC).

Scaling with the system size

The compute time scales as [math]\displaystyle{ N^4 }[/math] with the system size. The Dyson equation must be inverted at every frequency point (NOMEGA), so the total cost grows linearly with the number of frequency points.

Additional resources

Tutorials

How to

References

Pages in category "Time-dependent density functional theory"

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