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

Approximations to the interaction between electrons and holes are controlled in the INCAR by the tags LHARTREE, LADDER, and LFXC, which can be set to either .TRUE. or .FALSE.. Below we provide an explanation of what interaction term each tag controls.

{{{2}}}

Independent-particle approximation

In this approximation all interaction terms in the hamiltonian are turned off. You can do this by setting all tags (LHARTREE, LADDER, LFXC) to .FALSE.. This means that the computed spectrum will be equal to what you obtain from a ground-state calculation with LOPTICS=.TRUE.. However, you can use this calculation to test if everything is in order with your input files and the workflow is properly setup.

Hartree exchange potential

With the tag LHARTREE=.TRUE. the interaction terms in the hamiltonian include the unscreened Coulomb exchange. These terms are also known as the bubble diagrams from many-body perturbation theory (MBPT). If both LFXC and LADDER are set to .FALSE., this will be equivalent to running random phase approximation (RPA) calculation.

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.

The missing interaction between electrons and holes from either LFXC or LADDER means that you cannot properly describe bound excitons, which is a known problem of RPA. However, you can still use this setup to compute the electron energy-loss function, EELS. TODO: ADD LINK TO EELS IN CATEGORY:DIELECTRIC PROPERTIES ONCE IT'S ADDED TO THAT PAGE.

Screened two-particle interaction

Exchange-correlation effects from time-dependent density functional theory

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 accuracy of local exchange-correlation kernels when compared to experimental data 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 fail to properly reproduce the binding energies of electron-hole pairs in periodic systems or very large molecules.

Ladder diagrams from many-body perturbation theory

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

At the present, the screened interaction has to be computed from a model dielectric function.

LHFCALC = .TRUE.
LMODELHF = .TRUE.
HFSCREEN = 1.26
AEXX = 0.1

and it is also advisable to check the recommended values of HFSCREEN and AEXX for specific materials on the page of each respective tag.

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

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 the propagation is controlled by the NBANDSO and NBANDSV tags, 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. Due to the properties inherent to the Fourier transform used to integrate the time-dependent dipole moments, $T_\mathrm{max} \approx 1/\mathrm{CSHIFT}$ and $\Delta t \approx 1/\mathrm{OMEGAMAX}$.

The tag CSHIFT also controls 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}$.

Setting CSHIFT = 0.1 ~ 0.01 is often a good choice, as lower values will lead to unnecessarily long propagation times and spectra with very narrow peaks. OMEGAMAX is automatically from the maximum energy difference between occupied and unoccupied states, but you can lower this number to decrease the number of pairs used in the basis set.

Choosing the direction of perturbation

The dipole momentum vector $\mu^i$ is direction dependent, which means that you can choose along which direction the propagation takes place. In the INCAR this is controlled by IEPSILON, which can take values of 1, 2, 3 (corresponding to x, y, and z direction, respectively), and 4 (corresponding to all directions).

While choosing a single direction of propagation decreases the computing time, it is important that you pay attention to the symmetries of the material you are studying. For example, in the case of bulk silicon, since the material has cubic symmetry you only need to propagate along one direction (x, or y, or z). However, for a material like hexagonal boron nitride, the crystal symmetries make the x and y directions not equivalent. For this system you should propagate along both x and y and then take the average of the dielectric function along both directions to compute the absorption spectrum.