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.

Warning: All features presented in this tag are only available from VASP.6 or higher!

Macroscopic polarizability as a time-dependent integral

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

$\epsilon_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+c . c .\right) e^{-i(\omega-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 \hbar \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$.

The level of approximation is controlled by the Hamiltonian chosen in the time-dependent Schrödinger equation. 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.

By projecting the initial vector $\mu^i$ on $\xi^j(t)$ at each instant, VASP can compute the dielectric response at different directions ($i,j=x,y,z$). 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, for instance, while for cubic systems it is enough to compute $\epsilon_{xx}(\omega)$, for 2D systems such as hexagonal boron nitride it is advisable to compute both $\epsilon_{xx}(\omega)$ and $\epsilon_{yy}(\omega)$ and then analyse their average.

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 of the dipole-moment vector $\xi^j(t)$. Control over which valence and conduction states are added in time-dependent run is done with the NBANDSO and NBANDSV flags, respectively.

Users should note that while the inclusion of conduction states and the right potential makes comparison time-dependent results with data resulting from TDDFT or BSE calculations, 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.

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

By default, the time propagation code includes the Hartree and local-field effects (LHARTREE=.TRUE. and LFXC=.TRUE.). Results in the independent particle approximation can be calculated by setting both variables to .FALSE.. Please note that setting one variable to .FALSE. and the other to .TRUE. is not currently supported.

LHARTREE: activation of local field effects

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. 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 TDDFT 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 variables

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

Choice of bands variables

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$, while also automatically setting up the necessary number of steps to ensure numerical stability of the time-evolution 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

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

TIMEEV versus BSE and Casida TDDFT

It is important to distinguish the difference between time-evolution and what is done for the dielectric function coming from BSE and TDDFT calculations. The starting point comes from the fact that they are intended to describe different things. TIMEEV describes the evolution of the single-particle states and does not include interactions between electrons and holes. On the contrary, both BSE and TDDFT deal with the equation of motion of two-particle states and are, in principle, able to describe optical phenomena.

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 Tags and Sections

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

---