Time-evolution algorithm

The macroscopic dielectric function, $\epsilon_{ij}(\omega)$, measures how a given dielectric medium reacts when subject to an external electric field. From $\epsilon_{ij}(\omega)$ one can extract several optical properties such as absorption, optical conductivity, reflectance. However, it is important that the interacting electrons and holes are taken into account. This makes the evaluation of the macroscopic dielectric function more involved, since it goes beyond the single-particle level, working either at the Bethe-Salpeter or time-dependent density-functional theory level.

Within VASP, users can select two different methods for how $\epsilon_{ij}(\omega)$ is computed. The first is based on the eigendecomposition of the electron-hole Hamiltonian, $H^\mathrm{exc}$. It allows for the evaluation of $\epsilon_{ij}(\omega)$ by first obtaining the eigenvalues and eigenvectors of $H^\mathrm{exc}$ and it is based on the Bethe-Salpeter equation or the Casida equation. The second method transforms the mathematical expression of $\epsilon_{ij}(\omega)$ into a time-dependent integral. By propagating in time the dipolar moments and then applying a Fourier transform, it can compute $\epsilon_{ij}(\omega)$.

The advantage of the later method in comparison to the former is related to the cost, with the time-dependent integral being $O(N^2)$, while the eigendecomposition has a cost of $O(N^3)$, where $N$ is the rank of $H^\mathrm{exc}$. This means that for very large numbers of bands or k-points, the time-dependent formalism is cheaper than the eigendecomposition method.

Below we present a brief description of the method, from its theoretical support to how calculations should be performed, with the relevant approximations needed in the two-particle Hamiltonian.

The macroscopic-dielectric function as a time-dependent integral

The starting point is that one can re-write $\epsilon_{ij}(\omega)$ as a time-dependent integral. It starts from its expression, given by

$\epsilon^M(\omega)=1+\frac{4 \pi}{\Omega_0} \sum_\lambda\left|\sum_{c v \mathbf{k}} \mu_{c v \mathbf{k}} A_{c v \mathbf{k}}^\lambda\right|^2\left[\frac{1}{\omega+E_\lambda+\mathrm{i} \eta}-\frac{1}{\omega-E_\lambda+\mathrm{i} \eta}\right]$,

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 the dipolar moment associated to the the conduction $c$ and valence band $v$, and k-point $k$, $\lambda$ is the index of the eigenstate of $H^\mathrm{exc}$, with $A^\lambda$ and $E_\lambda$ being the associated eigenvector and eigenvalue. This equation can be brought into operational form,

$\epsilon^M(\omega)=1+\frac{4 \pi}{\Omega_0}\left\langle\mu\left|\left[\frac{1}{\omega+\mathrm{i} \eta+\hat{H}^{\mathrm{exc}}}-\frac{1}{\omega+\mathrm{i} \eta-\hat{H}^{\mathrm{exc}}}\right]\right| \mu\right\rangle$

by using the spectral decomposition $\left[\hat{H}^{\mathrm{exc}}-\omega\right]^{-1}=\sum_\lambda \frac{\left|A_\lambda\right\rangle\left\langle A_\lambda\right|}{E_\lambda-\omega}$. Then, one can bring the new expression of $\epsilon(\omega)$ into a time-dependent integral, by using

$\frac{1}{\omega+\mathrm{i} \eta-\hat{H}^{\mathrm{exc}}}|\mu\rangle=-\mathrm{i} \int_0^{\infty} e^{-\mathrm{i}\left(\omega-\hat{H}^{\mathrm{exc}}+\mathrm{i} \eta\right) t}|\mu\rangle=-\mathrm{i} \int_0^{\infty} e^{-\mathrm{i}(\omega+\mathrm{i} \eta) t} e^{\mathrm{i} \hat{H}^{\mathrm{exc}}t}|\mu\rangle$,

and recognising that $e^{\mathrm{i} \hat{H}^{\mathrm{exc}}t}|\mu\rangle = |\xi(t)\rangle$ is the exponential form of a time-dependent equation. These considerations allow the expression of $\epsilon_{ij}(\omega)$ to be written as

$\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 \eta) t}$,

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}^\mathrm{exc}(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$.

To compute the dielectric function with this method, VASP evaluates and stores at each time step the projections of $\left.\mid \xi^j(t)\right\rangle$ over $\left.\mid \mu^i\right\rangle$, $c^{ij}_{cv\mathbf k}(t) = \langle \mu^i_{cv\mathbf k}|\xi^i_{cv\mathbf k}(t)\rangle$. It is the fact that all these operations are of the matrix-vector type that makes this method having a cost of the order of $O(N^2)$.

Perturbing all transitions with a delta-like potential

explain how the system is perturbed with the delta-potential

The many-body terms in the hamiltonian

Screened two-particle interaction

Ladder diagrams from many-body perturbation theory

explain what different components are included in H (LFXC,LHARTREE,LADDER, or none)

Step-by-step instructions

Step 1: ground state with extra empty states

The starting point is a ground-state calculation which includes extra empty states, whose number is controlled in the INCAR file with the tag NBANDS. Taking the following input file as an example:

SYSTEM = Si
NBANDS = 12
ISMEAR = 0 ; SIGMA = 0.05
ALGO = N
LOPTICS = .TRUE.
KPAR = 4

8 empty bands are chosen (silicon has 4 occupied bands in the pseudo-potential file, thus making a total of 12 bands). However, with ALGO=N, VASP will employ an iterative diagonalization algorithm, meaning that the last conduction states will not be converged with the same accuracy as the occupied states. It is possible to avoid this by setting ALGO=Exact, or by increasing the number of bands to make sure that the states which will be used in the time-propagation step are converged with the same level of accuracy.

Finally, with LOPTICS=.TRUE., VASP will compute the dipole momentum for each possible $v\to c$ transition (recall the definition of the dipole momentum vector). These are written in the file WAVEDER, which will be used in the next step.

Step 2: time-evolution run

One the ground state with extra empty states is computed, you can copy the resulting WAVECAR and WAVEDER files to another directory where you will perform the time propagation. The following will be used as an example INCAR file:

SYSTEM = Si
ALGO = TIMEEV
!Information about the bands
NBANDS = 12
NBANDSO = 4
NBANDSV = 8
!Smearing parameters
ISMEAR = 0 ; SIGMA = 0.05
!Direction of propagation
IEPSILON = 1 
!Parallelization options
KPAR = 4
!Time-propagation parameters
NELMGW = 2000
CSHIFT = 0.1
OMEGAMAX = 20

We begin by noticing that now ALGO is set to TIMEEV, meaning that VASP will now perform a time-propagation calculation.

Setting up the bands

With NBANDS=12 you inform VASP that there are 12 states in total in the WAVECAR and WAVEDER. This must be consistent with Step 1! The number of occupied and unoccupied states that are used in th propagation is controlled by the NBANDSO and NBANDSV, respectively. To choose which bands to use, it is advisable to understand the type of property you want to compute. For instance, in the case of optical absorption, materials are probed within a few hundreds of milli-electronvolt of the band gap. In this case it means that you do not need to use states that lie too far away in energy from the band extrema.

In this example VASP will to use NBANDSO=4 occupied and NBANDSV=8 unoccupied states during the time propagation. You should be aware that there is no need to use the total number of bands set up by NBANDS, but still NBANDSO+NBANDSV cannot be larger than NBANDS.

Setting up the time-step

Since VASP is now integrating a time-dependent differential equation, you need to specify the time-step used to propagate the dipole moments. By default, VASP will use 20000 steps, however you can set a different number with the tag NELMGW. Nevertheless, NELMGW must be larger than 100, otherwise VASP will revert to the default.

This means that it is the length of the time step and the maximum propagation time that are dependent on the system in case and your input tags. These are controlled via the CSHIFT and OMEGAMAX tags.

The next parameter, CSHIFT controls the $\eta$ parameter, which is also the width used in the plotting of the dielectric function, since

$\frac{1}{\omega - E_\lambda + \mathrm i \eta} = \frac{1}{\omega - E_\lambda} - \mathrm i\pi \delta(\omega - E_\lambda)$

and the $\delta$-function is approximated as $\delta(\omega - E_\lambda) = \lim_{\eta\to 0^+}\frac{1}{\pi}\frac{\eta}{(\omega - E_\lambda)^2+\eta^2}$.