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)$ it is possible to 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 two-particle level via either the Bethe-Salpeter equation or time-dependent density-functional theory.

For both frameworks, BSE and TDDFT, users can select two different strategies to compute $\epsilon_{ij}(\omega)$. 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 initially obtaining the eigenvalues and eigenvectors of $H^\mathrm{exc}$ and then using both to evaluate $\epsilon_{ij}(\omega)$. This strategy is based on the Bethe-Salpeter equation or the Casida equation. The second strategy transforms the mathematical expression of $\epsilon_{ij}(\omega)$ into a time-dependent integral. By propagating the dipolar moments in time and then applying a Fourier transform, it can compute $\epsilon_{ij}(\omega)$.

The advantage of the later method in relation to the former is related to their cost. The time-dependent integral has a cost of the order $O(N^2)$, while the eigendecomposition has a cost of the order $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$, 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 definition of $\epsilon_{ij}(\omega)$ can be brought into operational form with

$\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}$.

The new expression of $\epsilon(\omega)$ can be brought into a time-dependent integral, using the fact that

$\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 recognizing that $e^{\mathrm{i} \hat{H}^{\mathrm{exc}}t}|\mu\rangle = |\xi(t)\rangle$ is the exponential form of a time-dependent equation operator. These considerations allow for 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 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

In order to probe all possible $v\to c$ transitions, a time-dependent term is added to the hamiltonian

$V_\mathrm{ext}(\mathbf r, t) = \lambda \mathbf r\cdot \mathbf D\delta(t),$

where $\lambda$ is the perturbation strength parameter and $\mathbf D$ is the electric displacement field. The narrow (in time) potential allows all bands in the occupied and unoccupied manifolds to be included in the transition space. The constant displacement field replicates the long wavelength limit (i.e. $\mathbf q \to 0$).

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.

Mind: The default setup for VASP is LHARTREE and LADDER set to .FALSE., while LFXC is set to .TRUE.. This means that if no tags are set in the INCAR the time-propagation run will using the TD-DFT kernel.

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.

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 from MBPT. The interaction hamiltonian will include the screened exchange interaction potential, $W(\mathbf r,\mathbf r';\omega)$, which has the correct long-range behaviour. This means that you can obtain bound electron-hole pairs in bulk materials and large molecules.

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

${\varepsilon}_{\mathbf{G},\mathbf{G}}^{-1}(\mathbf{q})=1-(1-{{\varepsilon}_{\infty}^{-1}})\text{exp}\left(-\frac{|\mathbf{q+G}|^2}{4{\lambda}^2}\right)$.

To activate this you must set both tags, LHFCALC and LMODELHF to .TRUE.. Also, you must provide VASP with both HFSCREEN ($\lambda$) and AEXX (${\varepsilon}_{\infty}^{-1}$) to control the range separation parameters in the model dielectric function.

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
!Particle interactions
LHARTREE = .TRUE.
LADDER = .TRUE.
LFXC = .FALSE.
LHFCALC = .TRUE.
LMODELHF = .TRUE.
AEXX = 0.088
HFSCREEN = 1.26

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 destroy the equivalency between the x and y directions. 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.

Comparison to other methods

VASP offers two other methods with which you can compute the macroscopic dielectric function. These are based on eigen-decomposition of the two particle hamiltonian, $H^\mathrm{exc}$. While more expensive than time-evolution, both these methods are able to compute eigenvalues and eigenstates of $H^\mathrm{exc}$, meaning that you can have direct access to the excitation energies of a system.

Bethe-Salpeter equation

Here the full Bethe-Salpeter equation is employed. The interaction hamiltonian is built using the dielectric function from RPA, and has the right behaviour in the long range regime. This means that it can accurately describe bound excitons in solids and large molecules. However, it is more costly than time-evolution, scaling with $N_\mathrm{rank}^3$.

Casida equation

Similar to the Bethe-Salpeter equation, the Casida equation employs an eigensolver method to compute the dielectric function. The key differences are that the Casida method does not require a preceding GW run to compute the RPA screening and can be performed with either DFT or hybrid-functional orbitals and energies.

Related tags and articles

NBANDSO

NBANDSV

IEPSILON

LHARTREE

LADDER

LFXC

LHFCALC

LMODELHF

AEXX

HFSCREEN

Time-dependent density-functional theory calculations

Bethe-Salpeter-equations calculations