Time-evolution algorithm

The macroscopic dielectric function, $\epsilon_{ij}(\omega)$, measures how a given dielectric medium reacts when subject to an external electric field. While VASP possesses several methods to compute this function, here we focus on a $O(N^2)$ strategy, which can be useful when dealing with large systems. This method comes from the fact that the computation of the dielectric function can be reformulated as a time-propagation problem, instead of an eigenvalue equation. Below we present a brief description of the method, as well as the relevant INCAR variables.

Mind: This feature presented here are only available for VASP.6 or higher!

Macroscopic polarizability as a time-dependent integral

The starting point is that one can re-write $\epsilon_{ij}(\omega)$ as a time-dependent integral

$\epsilon_{ij}(\omega)=\delta_{ij}-\frac{4\pi e^2}{\Omega}\int_0^{\infty} \mathrm{d} t \sum_{c,v,\mathbf{k}}\left(\langle\mu^j_{cv\mathbf{k}}| \xi^i_{cv\mathbf{k}}(t)\rangle+ \mathrm{c.c.}\right) e^{-\mathrm i(\omega-\mathrm i \delta) t}$,

where $\mu_{v c \mathbf{k}}^j=\frac{\left\langle c \mathbf{k}\left|v_j\right| v \mathbf{k}\right\rangle}{\varepsilon_c(\mathbf{k})-\varepsilon_v(\mathbf{k})}$ is a vector of all dipole moments. The fundamental aspect behind this transformation is that the new, time-dependent vector $\left.\mid \xi^j(t)\right\rangle$ follows the equation

$\mathrm i \frac{\mathrm d}{\mathrm d t}\left|\xi^j(t)\right\rangle=\hat{H}(t)\left|\xi^j(t)\right\rangle,$

with the initial vector elements given by $\left|\xi^j(0)\right\rangle=\left|\mu^j\right\rangle$. It is worthy to note that $H(t)$ is split into two parts, $H(t) = H_0 + V(t)$, where $H_0$ is the ground-state Hamiltonian, and $V(t)$ is the time-dependent perturbation.

To probe the dielectric response, $V(t)$ is set as proportional to a delta-pulse, i.e. $V(t) = \lambda\mathbf r\cdot\mathbf D\delta(t)$. With this, all possible transitions between the conduction and valence states included in the INCAR are probed during time-propagation. Control over which valence and conduction states are added in time-dependent run is done with the NBANDSO and NBANDSV flags, respectively.

By propagating the projection of $\mu^i$ onto $\xi^j(t)$ at each instant, $c_{cv\mathbf k}^{ij}(t)$, VASP can compute the dielectric response along different directions ($i,j=x,y,z$) by evaluating the Fourier transform of each time-dependent coefficient.

Mind: Since the problem is now a time-dependent equation, there is no computation of eigenvalues and eigenvectors of the two-particle Hamiltonian, unlike what happens in the Bethe-Salpeter formulation.

With the variable IEPSILON users can control the Cartesian direction along which the Dirac delta pulse is applied. The default value, IEPSILON = 4, performs three independent calculations for an electric field in x, y, and z direction, and thus is the most expensive. The choice of direction is system dependent: while for cubic systems it is enough to compute $\epsilon_{xx}(\omega)$, for 2D systems (eg. monolayer hexagonal boron nitride) it is advisable to compute both $\epsilon_{xx}(\omega)$ and $\epsilon_{yy}(\omega)$ and then analyse their average. This assuming that the $(001)$ direction is perpendicular to the plane in which the atoms are laid.

Users should note that a sufficient large number of such states must be included in a preceding ground-state calculation, where all the energy levels and dipole moments are computed. Note, however, that unoccupied orbitals are not propagated, which saves computational time. Also, if the interest is to compare time-evolution results with data coming from BSE/TDDFT calculations, users should set the same values for NBANDSO and NBANDSV in both calculations.

Finally, to perform this calculation users must set ALGO=TIMEEV in the INCAR file. A small example for bulk silicon is provided in the last section of this page.

Levels of approximation for the exchange-correlation potential

The interaction between electrons and holes can be controlled at three levels of approximation, depending on the input provided by the user. They are controlled by the INCAR variables, LADDER, LHARTREE, and LFXC. In this section we provide a brief description and some advice on when to use them.

LFXC: activation of the exchange-correlation kernel

This is the default setting if none of the aforementioned variables are included in the INCAR. Setting LFXC=.TRUE. includes the local exchange-correlation kernel, $f_\mathrm{xc}$ in the time-propagation

$f_{\mathrm{xc}}^{\text {loc }}\left(\mathbf{r}, \mathbf{r}^{\prime}\right)=\frac{\delta^2\left\{E_{\mathrm{c}}^{\mathrm{DFT}}+\left(1-c_{\mathrm{x}}\right) E_{\mathrm{x}}^{\mathrm{DFT}}\right\}}{\delta \rho(\mathbf{r}) \delta \rho\left(\mathbf{r}^{\prime}\right)},$

where $c_X$ controls the fraction of the exchange energy functional that is included in the kernel (see AEXX). This allows users to perform time-dependent calculations using hybrid functionals.

The use of local exchange-correlation kernels if often limited to molecules up to around 100-200 atoms. As these kernels often lack the long range component that goes as $-1/q^2$ (where $q$ is the momentum difference between the electron and the hole), they are can provide wrong results for periodic systems or very large molecules.

LADDER: activation of the screened exchange interaction

By setting LADDER=.TRUE., users can turn on the inclusion of the ladder diagrams, which account for what are usually called excitonic effects. These come from the screened exchange interaction potential, $W(\mathbf r,\mathbf r';\omega)$ and have the right long-range behaviour, allowing for the formation of bound electron-hole pairs in bulk materials.

For small molecules the inclusion of ladder diagrams is often unnecessary, but the long range interaction becomes important for large molecules.

LHARTREE: activation of the unscreened Coulomb interaction

With the flag LHARTREE=.TRUE. users can activate the inclusion of direct interaction terms in the Hamiltonian, also known as the bubble diagrams from many-body perturbation theory (MBPT). This is the equivalent of running a BSE calculation in the random phase approximation. Note that at the end, the dielectric function reported in the output files is the macroscopic dielectric function, where no contributions from local fields (i.e. terms with finite $\mathbf G$) are included.

LFXC: activation of the exchange-correlation kernel

Setting LFXC=.TRUE. includes the local exchange-correlation kernel, $f_\mathrm{xc}$ in the time-propagation

$f_{\mathrm{xc}}^{\text {loc }}\left(\mathbf{r}, \mathbf{r}^{\prime}\right)=\frac{\delta^2\left\{E_{\mathrm{c}}^{\mathrm{DFT}}+\left(1-c_{\mathrm{x}}\right) E_{\mathrm{x}}^{\mathrm{DFT}}\right\}}{\delta \rho(\mathbf{r}) \delta \rho\left(\mathbf{r}^{\prime}\right)},$

where $c_X$ controls the fraction of the exchange energy functional that is included in the kernel (see AEXX). This allows users to perform time-dependent calculations using hybrid functionals.

Choice of time-step

By default, the number of time steps is chosen automatically by VASP and and set to 20000. Alternatively, the number of time steps can be set directly by the tag NELM. In this case, the user-defined number of steps needs to be large than 100. Otherwise, the value of NELM will be discarded and VASP will revert to the default value.

Summary of required INCAR flags

Run-level

ALGO = TIMEEV - sets the calculation type to a time-evolution run.

Choice of bands

NBANDSO - number of occupied (valence) bands to be included.

NBANDSV - number of virtual (conduction) bands to be included.

Approximations to the potential

LHARTREE - include or exclude the Hartree exchange terms (i.e. bubble diagrams) in the Hamiltonian.

LFXC - include or exclude the $f_\mathrm{xc}$ kernel in the time-propagation.

AEXX - controls the fraction of exchange functional that is included in the $f_\mathrm{xc}$.

Time-propagation and spectrum

CSHIFT - controls the broadening of the peaks in the spectra of $\epsilon_{ij}(\omega)$, while also automatically setting up the maximum time used in the time-propagation algorithm.

IEPSILON - sets up the directions along which the dielectric function is computed.

NELM - sets the number of steps to be used in time-propagation. Users must set this number to be greater than 100 if they want to avoid the default value (20000).

OMEGAMAX - maximum frequency (in eV) for which the dielectric function is computed.

TIMEEV versus TDHF

VASP provides another algorithm with which users can compute the macroscopic dielectric function, TIMEEV. For TDHF, however, the problem is formulated in terms of an eigenvalue decomposition via the Casida equation. Both approaches are equivalent, although TIMEEV is faster for adiabatic LDA and PBE, while TDHF is faster for hybrid functionals.

Example

A typical calculation requires two steps. First, a ground state calculation:

SYSTEM = Si
NBANDS = 12 ! even 8 bands suffice for Si
ISMEAR = 0 ; SIGMA = 0.05
ALGO = N
LOPTICS = .TRUE.
KPAR = 4     ! assuming we run on 4 cores, this will be the fastest

Second, the actual time propagation:

SYSTEM = Si
NBANDS = 12  ! even 8 bands suffice for Si
ISMEAR = 0 ; SIGMA = 0.05
ALGO = TIMEEV
IEPSILON = 1 ! cubic system, so response in x direction suffices
NBANDSO = 4  ; NBANDSV = 8  ;  CSHIFT = 0.1
KPAR = 4     ! assuming we run on 4 cores, this will be the fastest

In this case, OMEGAMAX is set automatically to the maximal transition energy (about 25 eV in this example). Reducing the number of considered transitions, and thus reducing OMEGAMAX will increase both the duration of time steps and their number.

For standard DFT calculations, the time-propagation code is so fast that very dense k-point grids can often be used.

Related variables and articles

ALGO, CSHIFT, LHARTREE, LFXC, NBANDSV, NBANDSO, OMEGAMAX

---