GW and dielectric matrix

For GW calculations the frequency dependent dielectric matrix in the Random Phase Approximation (RPA) is determined via the polarizability ${\bf \chi}$ and the Coulomb potential $V$ using $\epsilon(\omega)= 1 -V \cdot \chi(\omega)$ .

Mind: low-scaling GW algorithms determine the dielectric matrix on the imaginary frequency axis and cannot be used to calculate

{\bf \epsilon}

on the real frequency axis.

The real-frequency dependent dielectric matrix can be calculated with the quartic-scaling GW implementation, which usability is limited to relatively small unit cells containing a few dozen atoms at maximum.

For instance, if one studies only the frequency dependent dielectric matrix, ALGO=CHI can be used and the calculation of GW quasi-particle energies is skipped. Here, two algorithms are selectable via LSPECTRAL and are discussed in the following.

Mind: A preceding DFT calculation with many unoccupied states is required for GW calculations.

Spectral Method

The spectral method is selected by default or manually with LSPECTRAL=.TRUE.. Here, VASP computes the polarizability in two steps. First, the spectral density of the polarizability is calculated with Fermi's golden rule^[1]^[2]

${\displaystyle S_{{\bf G G}'}({\bf q},\omega) = \sum_{nm}\sum_{\bf k q} \delta(\epsilon_{n\bf k} - \epsilon_{m\bf k-q} -\omega) \left\langle {n\bf k} \right| {\bf G} \left| {m\bf k-q} \right\rangle \left\langle {m\bf k-q}\right| {\bf G}' \left| {n\bf k}\right\rangle }$

The delta function is approximated with a triangular function centered at the transition energy and is rather efficient. The real and imaginary part of the polarizability is calculated in a second step using a Hilbert transformation^[2]

${\displaystyle \chi_{{\bf G G}'}({\bf q},\omega) = \int \mathrm{d}\omega' S_{{\bf G G}'}({\bf q},\omega') \left( \frac{1}{\omega-\omega' - i \eta } - \frac{1}{\omega+\omega' + i \eta } \right) }$

Here the value of η can be influenced by CSHIFT. However, this integration is performed semi-analytically by restricting the integration variable ω' to a frequency grid that is determined by NOMEGA, OMEGATL, OMEGAMIN and OMEGAMAX. Together with the approximation of the delta function in the spectral density (see above), the integration can be carried out analytically and one arrives essentially at a matrix-vector formula of the form

${\displaystyle \chi_{{\bf G G}'}({\bf q},\omega_j) = \sum_{k=1}^{\rm NOMEGA} t_{jk} S_{{\bf G G}'}({\bf q},\omega'_k) }$

Note, only the weights $t_{jk}$ depend essentially on the choice of CSHIFT. From the explicit form of these weights one can deduce that the smallest contributions to the dielectric function are terms with $\Delta_{\min}=\omega'_1-\omega'_2$ and are suppressed if $\eta>\Delta_\min$ .

Furthermore, the frequency grid $\omega_j$ of the complex polarizability is not restricted to the one of the integration variable. In fact, the resolution of the frequency grid of the polarizability can be set with NEDOS.

Tip: To reduce spurious contributions in the imaginary part of the dielectric function at small frequencies one can reduce CSHIFT and increase NOMEGA.

However, reducing CSHIFT in the spectral method introduces additional spikes in the dielectric function. Thus for visualization, one often performs a Lorentzian filter on the raw data, even though such a filter typically shifts peaks of the original dielectric function. There is no particular reason why the Lorentzian filter should be used; Gaussian smearing is also perfectly acceptable. Note, filtering is essentially cosmetics.

As an alternative approach, which already includes some sort of smoothing, one can calculate the dielectric function directly.

Direct calculation of the dielectric function

The polarizability can be calculated directly without using a Hilbert transformation via the formula of Adler and Wiser.^[3]^[4]

Here CSHIFT influences the peak width of the dielectric function directly.

This method is selected with LSPECTRAL=.FALSE. and requires less memory. However, the algorithm is also much slower that the spectral one, but might yield smoother dielectric functions than the spectral method.

Local field effects

For both method, local field effects in the dielectric function can be studied with LRPA=.TRUE..

Output

The imaginary and real part of frequency dependent dielectric function is always determined by the GW routine. It can be conveniently grepped from the file using the command (note two blanks between the two words)

grep " dielectric  constant" OUTCAR

The first value is the frequency (in eV) and the other two are the real and imaginary part of the trace of the dielectric matrix. Note that two sets can be found on the OUTCAR file.

After a successful calculation VASP writes following data into OUTCAR.

The head of the microscopic dielectric function (without local field effects): $\epsilon_{\rm mic}(\omega) = \lim_{{\bf q}\to 0} \epsilon_{\bf 00}({\bf q},\omega)$

Inverse macroscopic dielectric tensor: ${\displaystyle \frac{1}{ \hat {\bf q} \cdot \epsilon_\infty(\omega)\cdot \hat {\bf q} } = \lim_{{\bf q}\to 0} \left[\epsilon^{-1}\right]_{\bf 00}({\bf q},\omega)}$

The latter potentially includes local field effects depending on LRPA.

Note, OUTCAR contains only the frequency grid restricted to NOMEGA points. The complete frequency dependence (NEDOS frequency points) is written to vasprun.xml and will have following format:

 <dielectricfunction comment="HEAD OF MICROSCOPIC DIELECTRIC TENSOR (INDEPENDENT PARTICLE)">
  <imag>
   <array>
    <dimension dim="1">gridpoints</dimension>
    <field>energy</field>
    <field>xx</field>
    <field>yy</field>
    <field>zz</field>
    <field>xy</field>
    <field>yz</field>
    <field>zx</field>
    <set> 
     <r>     0.0000     0.0000     0.0000     0.0000     0.0000     0.0000     0.0000 </r>
     <r>     0.4627     0.0000     0.0000     0.0000     0.0000     0.0000     0.0000 </r>
     <r>     0.9250     0.0000     0.0000     0.0000     0.0000     0.0000     0.0000 </r>
     <r>     1.3866     0.0000     0.0000     0.0000     0.0000     0.0000     0.0000 </r>
     <r>     1.8472     0.0000     0.0000     0.0000     0.0000     0.0000     0.0000 </r>
     <r>     2.3065     0.0000     0.0000     0.0000     0.0000     0.0000     0.0000 </r>
        ...   ^          ^          ^          ^          ^          ^          ^   
              |          |          |          |          |          |          |
            energy     eps_xx     eps_yy     esp_zz     esp_xy     eps_yz     eps_xz

Comparison with k-p and density perturbation theory

The calculated microscopic frequency dependent dielectric function without local field effects corresponds to the same function obtained using LOPTICS=.TRUE., as well as the one obtained from density functional perturbation routines (LEPSILON=.TRUE.). In fact, it is guaranteed that the results are identical to those determined using a summation over conduction band states (LOPTICS). Differences for LSPECTRAL=.FALSE. must be negligible, and can be solely related to a different complex shift CSHIFT (defaults for CSHIFT are different in both routines). Setting CSHIFT manually in the INCAR file will remedy this issue. If differences prevail, it might be required to increase NEDOS. For LSPECTRAL=.TRUE. differences can arise, because (i) the GW routine uses less frequency points and different frequency grids than the optics routine or again (ii) from a different complex shift. Increasing NOMEGA should remove all discrepancies.

Mind: The GW routine is the only routine capable to include local field effects for the frequency dependent dielectric function.

Technical tips

If full GW calculations are not required, it is possible to greatly accelerate the calculations, by calculating the response functions only at the $\Gamma$ -point. This can be achieved by setting the following values in the INCAR file:

 NKREDX = number of k-points in direction of first lattice vector
 NKREDY = number of k-points in direction of second lattice vector
 NKREDZ = number of k-points in direction of third lattice vector

The calculation of the QP shifts can be bypassed by setting ALGO=CHI. Furthermore, if only the static response function is required the number of frequency points should be set to NOMEGA=1 and LSPECTRAL=.FALSE.

[gajdos:prb:2006-1] M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, Phys. Rev. B 73, 045112 (2006).

[shishkin:prb:2006-2] M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006).

[adler:1962-3] S. L. Adler, Phys. Rev. 126, 413 (1962)

[wiser:1963-4] N. Wiser, Phys. Rev. 129, 62 (1963)

[1]

[2]

[3]

[4]