Time-dependent density-functional theory calculations

VASP offers a powerful module for performing time-dependent density-functional theory (TDDFT) or time-dependent Hartree-Fock (TDHF) calculations by solving the Casida equation . This approach can be used for obtaining the frequency-dependent dielectric function with the excitonic effects and can be based on the ground-state electronic structure in the DFT, hybrid-functional or even GW approximations. You can watch a lecture covering TDDFT theory and calculations on our YouTube channel.

From TDDFT (or TDHF) to Casida's equation

In time-dependent density functional theory the electrons follow a time-dependent Schrödinger equation

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

where all quantities are now dependent on a time-evolving density, [math]\displaystyle{ n(t) }[/math]. The external potential [math]\displaystyle{ V_\mathrm{ext}(t) }[/math] now is explicitly dependent on time, and which is used to perturb the system away from the ground state. For time-dependent Hartree-Fock (TDHF) one would only need to replace the exchange correlation potential, [math]\displaystyle{ V_\mathrm{xc} }[/math], with the Fock exchange term, [math]\displaystyle{ V_\mathrm{x} }[/math].

From linear-response theory it is possible to evaluate the change in the ground-state density due to the external potential,

[math]\displaystyle{ \chi(\mathbf 1,2) = \frac{\delta n(1)}{\delta V_\mathrm{ext}(2)}, }[/math]

or due to the total Kohn-Sham potential

[math]\displaystyle{ \chi_0(\mathbf 1,2) = \frac{\delta n(1)}{\delta V_\mathrm{Ks}(2)}. }[/math]

Both quantities are related via a Dyson equation. It is possible to show that

[math]\displaystyle{ \frac{\delta V_\mathrm{KS}(1)}{\delta V_\mathrm{ext}(2)} = \delta(1,2) + \frac{\delta(t_1-t_2)}{|\mathbf r_1 - \mathbf r_2|} + f_\mathrm{xc}(1,2), }[/math]

with [math]\displaystyle{ f_\mathrm{xc}(1,2) = \frac{\delta V_\mathrm{xc}(1)}{\delta n(2)} }[/math] being the exchange-correlation kernel. This quantity is dependent on the choice taken to the exchange-correlation energy functional. The application of a chain rule to [math]\displaystyle{ \chi(1,2) }[/math] leads to a Dyson equation

[math]\displaystyle{ \chi(1,2) = \chi_0(1,2) + \chi_0(1,3)\left[\frac{\delta(t_3-t_4)}{|\mathbf r_3 - \mathbf r_4|} + f_\mathrm{xc}(3,4)\right]\chi(4,2) }[/math]

While relatively simple, it is easier to solve this equation working within the frequency domain and in a basis set that considers transitions from valence to conduction states at the same k-point. The latter choice comes from the fact that, when dealing with optical measurements (e.g. absorption), we deal with neutral excitations and the momentum of the absorbed photons is practically zero. Once all of this is taken into account, the equation for [math]\displaystyle{ \chi(1,2) }[/math] is recast as an eigenvalue problem

[math]\displaystyle{ \left( \begin{matrix} A & B \\ -B^* & -A^* \end{matrix} \right) \left( \begin{matrix} X \\ Y \end{matrix} \right) = \omega \left( \begin{matrix} X \\ Y \end{matrix} \right), }[/math]

which is also known as the Casida equation. The A and B matrices are given by

[math]\displaystyle{ A_{vc}^{v'c'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'} + \langle cv'|\frac{1}{|\mathbf r_1-\mathbf r_2|}|vc'\rangle - \langle cv'|f_\mathrm{xc}(1,2)|c'v\rangle, }[/math]

while the B matrix is computed via

[math]\displaystyle{ B_{vc}^{v'c'} = \langle vv'|\frac{1}{|\mathbf r_1-\mathbf r_2|}|cc'\rangle - \langle vv'|f_\mathrm{xc}(1,2)|c'c\rangle. }[/math]

The A matrix describes the resonant (excitations) and anti-resonant (de-excitation) transitions, while the B matrix deals with the coupling between both. Due to the presence of this coupling, Casida's matrix is non-Hermitian.

This formalism is very similar to the Bethe-Salpeter equation, only now instead of the exchange-correlation self-energy, [math]\displaystyle{ \Sigma_\mathrm{xc} }[/math], and the screened interaction, [math]\displaystyle{ W }[/math], it deals with the exchange-correlation potential, [math]\displaystyle{ V_\mathrm{xc} }[/math], and with the exchange-correlation kernen, [math]\displaystyle{ f_\mathrm{xc} }[/math].

Solving Casida equation

The algorithm for solving the Casida equation can be selected by setting ALGO = TDHF. This approach is very similar to BSE but differs in the way the screening of the Coulomb potential is approximated. The TDHF approach uses the exchange-correlation kernel [math]\displaystyle{ f_{\rm xc} }[/math], whereas BSE requires the [math]\displaystyle{ W(\omega \to 0) }[/math] from a preceding GW calculation. Thus, in order to perform a TDHF calculation, one only needs to provide the ground-state orbitals (WAVECAR) and the derivatives of the orbitals with respect to [math]\displaystyle{ k }[/math] (WAVEDER).

Mind: Unlike BSE, TDHF calculations do not require [math]\displaystyle{ W(\omega \to 0) }[/math], i.e., Wxxxx.tmp

In summary, both TDHF and BSE approaches require a preceding ground-state calculation, however, the TDHF does not need the preceding GW and can be performed with the DFT or hybrid-functional orbitals and energies.

Time-dependent Hartree-Fock

The TDHF calculations can be performed in two steps:

1. ground-state calculation
2. optical absorption calculation

For example, an optical absorption calculation of bulk Si can be performed using a dielectric-dependent hybrid-functional described in Refs..

SYSTEM    = Si
ISMEAR    = 0 
SIGMA     = 0.05
NBANDS    = 16      ! or any larger desired value
ALGO      = D       ! Damped algorithm often required for HF type calculations, ALGO = Normal might work as well
LHFCALC   = .TRUE. 
LMODELHF  = .TRUE. 
AEXX      = 0.083
HFSCREEN  = 1.22
LOPTICS   = .TRUE.  ! can also be done in an additional intermediate step

In the second step, the dielectric function is evaluated by solving the Casida equation

SYSTEM    = Si
ISMEAR    = 0 
SIGMA     = 0.05
NBANDS    = 16     
ALGO      = TDHF
IBSE      = 0
NBANDSO   = 4       ! number of occupied bands
NBANDSV   = 8       ! number of unoccupied bands
LHARTREE  = .TRUE.
LADDER    = .TRUE.
LFXC      = .TRUE.
LMODELHF  = .TRUE. 
AEXX      = 0.083
HFSCREEN  = 1.22

THDF calculations can be performed for non-spin-polarized, spin-polarized, and noncollinear cases, as well as the case with spin-orbit coupling. There is, however, one caveat. The local exchange-correlation kernel is approximated by the density-density part only. This makes predictions for spin-polarized systems less accurate than for non-spin-polarized systems.

Time-dependent DFT calculation

The TDDFT calculation using the PBE exchange-correlation kernel can be performed by disabling the ladder diagrams LADDER = .FALSE., i.e., only the PBE exchange-correlation kernel is present in the Hamiltonian.

SYSTEM    = Si
ISMEAR    = 0 
SIGMA     = 0.05
NBANDS    = 16     
ALGO      = TDHF
IBSE      = 0
NBANDSO   = 4       ! determines how many occupied bands are used
NBANDSV   = 8       ! determines how many unoccupied (virtual) bands are used
LFXC      = .TRUE.
LHARTREE  = .TRUE.
LADDER    = .FALSE.

Mind: In TDDFT calculation, where the ladder diagrams are not included (LADDER=.FALSE.) or the fraction of exact exchange in the kernel is zero (AEXX=0), the resulting dielectric function lacks the excitonic effects.

VASP tries to use sensible defaults, but it is highly recommended to check the OUTCAR file and make sure that the right bands are included. The tag OMEGAMAX specifies the maximum excitation energy of included electron-hole pairs and the pairs with the one-electron energy difference beyond this limit are not included in the Hamiltonian.

The calculated frequency-dependent dielectric function, transition energies and oscillator strength values are stored in the vasprun.xml file.

Calculations beyond Tamm-Dancoff approximation

Calculations beyond Tamm-Dancoff approximation can be performed in the same manner as in the BSE.

Calculations at finite wavevectors

Calculations at finite wavevectors can be performed in the same manner as in the BSE.

---