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^[1]

{\displaystyle \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

{\displaystyle \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^[2], i.e. $V(t) = \lambda\mathbf r\cdot\mathbf D\delta(t)$ ^[3]. 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

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.

LADDER: activation of the screened exchange interaction

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

{\displaystyle 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.

References

[schmidt:prb:2003-1] W. G. Schmidt, S. Glutsch, P. H. Hahn, and F. Bechstedt, Efficient O(N2) method to solve the Bethe-Salpeter equation, Phys. Rev. B 67, 085307 (2003)

[kubo:jpsj:1957-2] R. Kubo, Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems, J. Phys. Soc. Jpn. 12, 570 (1957).

[sander:jcp:2017-3] T. Sander, G. Kresse, Macroscopic dielectric function within time-dependent density functional theory—Real time evolution versus the Casida approach , J. Chem. Phys. 146, 064110 (2017)

[1]

[2]

[3]