3.42 Linear macroscopic dielectric function and Kubo-Greenwood transport

The linear macroscopic dielectric tensor $\epsilon_{ij}(\omega)$ , the ratio between the average of the total potential in one unit cell and the external field, within the RPA is calculated. For derivation from the microscopic dielectric function $\epsilon(\vec{r},t;\vec{r^{\prime}},t^{\prime})$ and references see [7] (Chapter 1, 2 and appendix). From the complex frequency dependent dieletric function all other optical constants can be determined, e.g. optical conductivity, Loss function, Reflectivity, etc. Also see Appendix K in Ashcroft/Mermin: solid state physics.
The frequency dependent real and imaginary part of the inter- (eq. 3.190) and intra- (eq. 3.189) band contribution to the linear dielectric tensor is calculated as post-processing after convergence of the SCF-cycle from (atomic units):

	$\displaystyle\epsilon_{ij}(\omega)=\delta_{i,j}$	$\displaystyle-\frac{4\pi}{V_{cell}\omega^{2}}\sum_{n,\mathbf{k}}\left(-\frac{% \partial f_{0}}{\partial\epsilon}\right)_{\epsilon_{n,\mathbf{k}}}p_{i;n,n,% \mathbf{k}}p_{j;n,n,\mathbf{k}}$		(3.189)
		$\displaystyle+\frac{4\pi}{V_{cell}}\sum_{\mathbf{k}}\sum_{c,v}\frac{p_{i;c,v,% \mathbf{k}}p_{j;c,v,\mathbf{k}}}{\left(\epsilon_{c,\mathbf{k}}-\epsilon_{v,% \mathbf{k}}-\omega\right)\left(\epsilon_{c,\mathbf{k}}-\epsilon_{v,\mathbf{k}}% \right)^{2}}$		(3.190)

$f$ is the Fermi-function and $V_{cell}$ the unit cell volume. The intra band part is singular at $\omega=0$ . Here the plasma frequency $\omega_{pl;ij}$ is calculated:

\omega^{2}_{pl;ij}=\frac{1}{\pi}\sum_{n}\int_{\vec{k}}p_{i;n,n,\vec{k}}p_{j;n,% n,\vec{k}}\delta(\epsilon_{n,\vec{k}}-\epsilon_{F})

(3.191)

NB:

The quantity $\omega^{2}_{\mathrm{pl};ij}$ produced by this output is known to be buggy and may produce incorrect values. Please refer to the correct implementation of the plasma frequency provided in the subsection below. For details see:

•

M. K. Pogodaeva, S. V. Levchenko, V. P. Drachev, and I. R. Gabitov, “Dielectric function of six elemental metals,” J. Phys.: Conf. Ser. 1890, 012008 (2021), https://doi.org/10.1088/1742-6596/1890/1/012008.IOP Publishing.

By adopting a Drude-like shape for the intra-band contributions a lifetime broadening $\Gamma$ is introduced and $\epsilon_{ij}^{inter}(\omega)$ becomes:

	$\displaystyle\mathrm{Im}(\epsilon_{ij}^{intra}(\omega))=\frac{\Gamma\omega^{2}% _{pl;ij}}{\omega(\omega^{2}+\Gamma^{2})}$		(3.192)
	$\displaystyle\mathrm{Re}(\epsilon_{ij}^{intra}(\omega))=1-\frac{\omega^{2}_{pl% ;ij}}{\omega^{2}+\Gamma^{2}}$		(3.193)

In the case of a spin unpolarized calculation eq. 3.190 and eq. 3.189 have to be multiplied by $2$ to yield the correct occupation.
The following quantities are needed/calculated:

-

$\epsilon_{n,\vec{k}}$ - the eigenvalue of the Kohn-Sham eigenstate $(n,\vec{k})$ .
-

$\delta(\epsilon_{n^{{}^{\prime}},\vec{k}}-\epsilon_{n,\vec{k}}-\omega)=\frac{1% }{\sqrt{2\pi}width}exp\left(-\frac{1}{2}\frac{(\epsilon_{n^{{}^{\prime}},\vec{% k}}-\epsilon_{n,\vec{k}}-\omega)^{2}}{\Gamma^{2}}\right)$ - Gauss function with width $\Gamma$ = $w i d t h$ for calculating the plasma frequency. In eq 3.190 $\omega$ is replaced by $\omega+\mathrm{i}\Gamma$ , introducing Lorentzian broadenig.
-

$p_{j;n^{{}^{\prime}},n,\vec{k}}=\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|% \nabla_{j}|\psi_{n\vec{k}}\right\rangle$ - the momentum matrix elements calculated form the real space basis functions in k-space and KS-eigenstate basis

\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|\nabla_{j}|\psi_{n\vec{k}}\right% \rangle=\sum_{ij}c_{in^{{}^{\prime}}}^{*\vec{k}}c_{jn}^{\vec{k}}\sum_{\vec{N},% \vec{M}}\mathrm{e}^{\mathrm{i}\vec{k}\left[\vec{T}\left(\vec{N}\right)-\vec{T}% \left(\vec{M}\right)\right]}\left\langle\phi_{i\vec{M}}|\nabla_{j}|\phi_{j\vec% {N}}\right\rangle

(3.194)

\left\langle\phi_{i\vec{M}}|\nabla|\phi_{i\vec{N}}\right\rangle=\int_{\text{% unit cell}}\mathrm{d}^{3}r\phi_{i,\vec{M}}\vec{\nabla}\phi_{j,\vec{N}}

(3.195)

with the real space functions $\phi_{i}(\vec{r})$ centered in unit cells shifted by $\vec{T}\left(\vec{N}\right),~{}\vec{N}=(N_{1},N_{2},N_{3})$ , see [36] for details.

Building up on the expressions for the dielectric function also the corresponding Kubo-Greenwood transport properties expressed via the Onsager coefficients $L_{ij}$ :

L_{ij}(\omega)=\frac{2\pi e^{4-i-j}}{3Vm^{2}_{e}\omega}\sum_{\bm{k},m,n}|% \langle\Psi_{m}|\bm{\hat{p}}|\Psi_{n}\rangle|^{2}\cdot\left(\frac{\epsilon_{% \bm{k}m}+\epsilon_{\bm{k}n}}{2}-\epsilon_{Fermi}\right)^{i+j-2}\cdot\left(f_{% \bm{k}m}-f_{\bm{k}n}\right)\delta\left(\epsilon_{\bm{k}n}-\epsilon_{\bm{k}m}-% \hbar\omega\right)

(3.196)

In this representation $L_{11}$ corresonds to the optical conductivity $\sigma(\omega)$ . Furthermore the (frequency dependent) Seebeck coefficient is easily obtainable:

S=\frac{L_{12}}{TL_{11}}

(3.197)

3.42.1 Plasma frequency

Definition

The plasma frequency tensor is computed as

\bar{\omega}_{ij}^{2}=\bar{\omega}_{ij}^{\uparrow 2}+\bar{\omega}_{ij}^{% \downarrow 2}=-\frac{4\pi}{V}\sum_{\sigma n\mathbf{k}}\frac{\partial f\left(% \varepsilon_{n}^{\sigma}\right)}{\partial\varepsilon_{n}^{\sigma}}\left(% \mathbf{e}_{i}\frac{\partial\varepsilon_{n}^{\sigma}(\mathbf{k})}{\partial% \mathbf{k}}\right)\left(\mathbf{e}_{j}\frac{\partial\varepsilon_{n}^{\sigma}(% \mathbf{k})}{\partial\mathbf{k}}\right),

(3.198)

where $f(\varepsilon_{n}^{\sigma})$ is the Fermi function for spin $\sigma$ and band $n$ , $V$ is the unit cell volume, $\varepsilon_{n}^{\sigma}(\mathbf{k})$ are Kohn–Sham eigenvalues, and $\mathbf{e}_{i}$ is a unit vector in Cartesian direction $i$ .

Kohn–Sham framework

In periodic boundary conditions, the Kohn–Sham equations are $\mathbf{k}$ -dependent:

\sum_{\beta}c_{\beta n}^{\sigma}(\mathbf{k})\sum_{\mathbf{R}}\left(h_{\alpha% \beta}^{\sigma\mathbf{R}}-\varepsilon_{n}^{\sigma}(\mathbf{k})s_{\alpha\beta}^% {\mathbf{R}}\right)e^{i\mathbf{k}\mathbf{R}}=0,

(3.199)

where $c_{\beta n}^{\sigma}(\mathbf{k})$ are k-dependent coefficients, and $h_{\alpha\beta}^{\sigma\mathbf{R}}$ , $s_{\alpha\beta}^{\mathbf{R}}$ are real-space Hamiltonian and overlap matrix elements.

The derivative of a Kohn–Sham eigenenergy is then

\frac{d\varepsilon_{n}^{\sigma}(\mathbf{k})}{d\mathbf{k}}=\sum_{\alpha\beta% \mathbf{R}}i\mathbf{R}\,c_{\alpha n}^{\sigma*}(\mathbf{k})c_{\beta n}^{\sigma}% (\mathbf{k})e^{i\mathbf{k}\mathbf{R}}\left[h_{\alpha\beta}^{\sigma\mathbf{R}}-% s_{\alpha\beta}^{\mathbf{R}}\varepsilon_{n}^{\sigma}(\mathbf{k})\right].

(3.200)

Spin–orbit coupling (SOC)

When the SOC Hamiltonian is included, spin-up and spin-down channels are mixed. The derivative becomes

\frac{d\varepsilon_{n}(\mathbf{k})}{d\mathbf{k}}=\sum_{\mathbf{R}}i\mathbf{R}e% ^{i\mathbf{k}\cdot\mathbf{R}}\sum_{\alpha,\beta}\left\{\sum_{\sigma}C_{\alpha n% }^{\sigma*}(\mathbf{k})C_{\beta n}^{\sigma}(\mathbf{k})\left[h_{\alpha\beta}^{% \sigma\mathbf{R}}-s_{\alpha\beta}^{\mathbf{R}}\varepsilon_{n}\right]+\sum_{% \sigma\sigma^{\prime}}C_{\alpha n}^{\sigma*}(\mathbf{k})C_{\beta n}^{\sigma^{% \prime}}(\mathbf{k})v_{\alpha\beta}^{\sigma\sigma^{\prime}\mathbf{R}}\right\},

(3.201)

where $v_{\alpha\beta}^{\sigma\sigma^{\prime}\mathbf{R}}$ are SOC matrix elements. Here, $C_{\alpha n}^{\sigma}(\mathbf{k})$ are SOC-perturbed eigenvectors expanded in the basis of unperturbed ones. In practice, if the SOC perturbation is small, terms involving $\partial C_{\alpha n}^{\sigma}(\mathbf{k})/\partial\mathbf{k}$ can be neglected.

Convergence requirements

Accurate plasma frequency calculations require dense k-point sampling. Typical guidelines:

•

Simple metals: $~{}150$ of k-points per direction for fcc Al
•

Noble metals (e.g. Au): more than 70 k-points per direction are typically required for convergence within 1–2%. For more details see:
•

M. K. Pogodaeva, S. V. Levchenko, and V. P. Drachev, “Spin-dependent plasma frequency from all-electron ab initio calculations including spin-orbit coupling,” Phys. Rev. B 107, 045113 (2023), https://doi.org/10.1103/PhysRevB.107.045113.
•

Complex systems like supercells with defects or metals with low symmetry may suffice < 10 k-points.
•

A. Tselin, M. Pogodaeva, S. Levchenko, A. Kalmykov, K. Garbuzov, A. Smirnov, and V. Drachev, “Exploring the dielectric function of platinum,” Phys. Rev. B 110, 195130 (2024), https://doi.org/10.1103/PhysRevB.110.195130.

Refer to caption — Figure 3.12: Example convergence behavior of the plasma frequency for gold with respect to the k-point sampling.

Tags for general section of control.in:

Tag: compute_dielectric(control.in)

Usage: compute_dielectric $\omega_{max}$ $n_{\omega}$
Purpose: Sets basic parameters for calculating the imaginary and real part of the inter-band and intra-band contribution to the linear dielectric tensor within the random-phase approximation. This keyword is specified once in control.in to set parameters for the dielectric tensor calculation.
$\omega_{max}$ is a real valued number. This defines the upper frequency limit (in eV) up to which which the dielectric function should be written. Choose 10 eV or larger in order to ensure that later transformations, such as Kramers-Kronig, are accurate.
$n_{\omega}$ is an integer number, i.e., the number of frequency points for which the dielectric function is calculated. This number should be at least 100 or, perhaps, much larger.

By setting this keyword, all dielectric tensor components (referenced to cartesian coordinates: $x x$ , $y y$ , $z z$ , $x y$ , $x z$ , $y z$ ) will be written to files. Additionally, the corresponding absorption coefficients (also in cartesian coordinates: $x x$ , $y y$ , $z z$ ) will be written as derived from the diagonal elements of the dielectric function. Note that these are really cartesian coordinates – they are not oriented relative to the lattice vectors in geometry.in.

The calculation is very sensitive to the k-point grid. An extremely high number of k-points may be needed for convergence of the dielectric function for systems with strong band dispersion or for metals. Please test the k_grid convergence carefully.

The momentum matrix elements of pairs of occupied and unoccupied states in the energy window [VBM-( $\omega_{max}$ + 10.0 eV), CBM+( $\omega_{max}$ + 10.0 eV)] (in eV) relative to the electronic chemical potential determined alongside the occupation numbers of the DFT single-particle states will be summed over.

The default broadening type and broadening width used for the delta function is 0.1 eV if Lorentzian broadening is used. However, due to a long and (in our view) spurious tail of the Lorentzian broadening into the band gap of semiconductors, we recommend to use Gaussian broadening instead. These settings can be changed via the dielectric_broadening keyword or they can be changed by just using the output dielectric keyword as shown in the example below.

In order to avoid numerical integration errors caused by including 0 eV energy, in the energy range, the minimum energy ( $\omega_{min}$ ) is instead automatically set to a specific value corresponding to the broadening width that was specified ( $\omega_{min}$ = broadening width / 10.0).

The resulting output files will be named dielectric_function_(directions).out (e.g., dielectric_function_x_x.out for the $x x$ tensor component) and absorption_(directions).out.

In order to test the impact of different broadening type and broadening parameters, the individual tensor component for the dielectric constants can also be specified via the output dielectric keyword. By setting this keyword, the code will only output the dielectric functions in the directions listed in the output dielectric keyword, instead of automatically writing all tensor components.

An example snippet for control.in that tests multiple dielectric function components and broadening methods in a single FHI-aims run is given below. Note that small broadening parameters require denser k-grids for numerical convergence:

compute_dielectric   10.0  1000

output dielectric  gaussian   0.1   x  x
output dielectric  gaussian   0.1   y  y
output dielectric  gaussian   0.1   z  z
output dielectric  gaussian   0.1   x  y
output dielectric  gaussian   0.1   x  z
output dielectric  gaussian   0.1   y  z

output dielectric  gaussian   0.07  x  x
output dielectric  gaussian   0.07  y  y
output dielectric  gaussian   0.07  z  z
output dielectric  gaussian   0.07  x  y
output dielectric  gaussian   0.07  x  z
output dielectric  gaussian   0.07  y  z

output dielectric  lorentzian 0.1   x  x
output dielectric  lorentzian 0.1   y  y
output dielectric  lorentzian 0.1   z  z
output dielectric  lorentzian 0.1   x  y
output dielectric  lorentzian 0.1   x  z
output dielectric  lorentzian 0.1   y  z

Tag: dielectric_broadening(control.in)

Usage: dielectric_broadening broadening_method width

Purpose: Changing the broadening function type and broadening parameters used in the dielectric calculation. To use this keyword, the compute_dielectric keyword must be specified in control.in
The Delta-distribution is approximated by a broadening function specified by the broadening_method option with a defined width specified by the width option (in eV). Gaussian (gaussian) and lorentz (lorentzian) broadenings are supported.

Tag: output dielectric(control.in)

Usage: output dielectric broadening_method width i j

Purpose: Output the ij tensor components (choices: i,j = x, y, z) of the imaginary and real parts of the inter-band and intra-band contribution to the linear dielectric tensor. This keyword may be specified multiple times in control.in to output more than one tensor component (possibly with different broadenings) per calculation. The RPA approximation (i.e. Lindhard theory) is used. To use this keyword, the compute_dielectric keyword must be specified in control.in.

Note: This keyword are just setted for testing purpose. By setting this keyword, only the directions you listed will be output. The Delta-distribution is approximated by a broadening function specified by the broadening_method option with a defined width specified by the width option (in eV). Gaussian (gaussian) and Lorentz (lorentzian) broadenings are supported. Good starting choices for parameters are broadening_method = gaussian and width= $0.05eV$ . We encourage the user to try out different broadenings by specifying this keyword multiple times.

Tag: compute_absorption(control.in)

Usage: compute_absorption width $E_{min}$ $E_{max}$ $\omega_{min}$ $\omega_{max}$ $n_{\omega}$ i use_gauss
Purpose: Calculate and output the i component (choices: x, y or z) of the linear absorption.

\alpha_{i}(\omega)=\frac{8\pi^{2}}{\omega V_{cell}}\sum_{c,v}\sum_{\vec{k}}% \left|p_{i;c,v,\vec{k}}\right|^{2}\delta(\epsilon_{c,\vec{k}}-\epsilon_{v,\vec% {k}}-\omega)\mathrm{d}\vec{k}

(3.202)

The momentum matrix elements in the energy window [ $E_{min}$ , $E_{max}$ ] (in eV) relative to the internal zero are summed up. The Delta-distribution is represented by a Gaussian function (use_gauss = .true.) or a Lorentz function (use_gauss = .false.) with width width (in eV) for $n_{\omega}$ $\omega$ -values in the interval [ $\omega_{min}$ , $\omega_{max}$ ] (in eV). The output file is named absorption_i.out. A good choice is: width= $0.1eV$ , usually it is enough to include only a few states below and above the fermi level in the energy window [ $E_{min}$ , $E_{max}$ ]. The unit is $a_{0}^{-1}$ . $V_{cell}$ the unit cell volume. (see page 38, eq. 3.13 and 3.14 of http://www.tddft.org/bmg/files/papers/85619.pdf, [39])
The calculation is very sensitive to the k-point grid, an extremely high number of k-points might be needed for convergence, especially for metals.

Tag: compute_momentummatrix(control.in)

Usage: compute_momentummatrix $E_{min}$ $E_{max}$ k-point
Purpose: Calculate and output the momentum matrix elements $\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|\nabla_{x}|\psi_{n\vec{k}}\right\rangle$ , $\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|\nabla_{y}|\psi_{n\vec{k}}\right\rangle$ , $\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|\nabla_{z}|\psi_{n\vec{k}}\right\rangle$ for the k-point with the number k-point for KS-eigenstates that are within the energy window [ $E_{min}$ , $E_{max}$ ] (in eV) relative to the internal zero. The output file is named element_k_k-point.dat. The unit is $a_{0}^{-1}$ (bohr^-1).
If you are conducting a cluster calculation (no periodic boundary conditions) make sure to set k-point to 1.
Setting k-point to 0 will output the momentum matrix elements for all k-points into one container file (mommat.h5). To use this functionality, FHI-aims has to be compiled with the external hdf5 module (see Sec. H.2 for details).

Tag: compute_dipolematrix(control.in)

Usage: compute_dipolematrix $E_{min}$ $E_{max}$ k-point
Purpose: Calculate and output the dipole matrix elements $\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|x|\psi_{n\vec{k}}\right\rangle$ , $\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|y|\psi_{n\vec{k}}\right\rangle$ , $\left\langle\psi_{n^{{}^{\prime}}\vec{k}}|z|\psi_{n\vec{k}}\right\rangle$ (position operator!) for the k-point with the number k-point for KS-eigenstates that are within the energy window [ $E_{min}$ , $E_{max}$ ] (in eV) relative to the internal zero. The output file is named element_k_k-point.dat. The unit is $a_{0}$ (bohr).
If you are conducting a cluster calculation (no periodic boundary conditions) make sure to set k-point to 1.
Setting k-point to 0 will output the dipole matrix elements for all k-points into one container file (dipmat.h5). To use this functionality, FHI-aims has to be compiled with the external hdf5 module (see Sec. H.2 for details).

Tag: compute_dipolematrix_k_k(control.in)

Usage: compute_dipolematrix_k_k $E_{min}$ $E_{max}$ k_k_method
Purpose: Calculate and output the dipole matrix elements $\left\langle\psi_{n^{{}^{\prime}}\vec{k^{\prime}}}|x|\psi_{n\vec{k}}\right\rangle$ , $\left\langle\psi_{n^{{}^{\prime}}\vec{k^{\prime}}}|y|\psi_{n\vec{k}}\right\rangle$ , $\left\langle\psi_{n^{{}^{\prime}}\vec{k^{\prime}}}|z|\psi_{n\vec{k}}\right\rangle$ (position operator!) for all (k, k’)-points for KS-eigenstates that are within the energy window [ $E_{min}$ , $E_{max}$ ] (in eV) relative to the internal zero. The output file is named dipmat_k.h5. The unit is $a_{0}$ (bohr). With the option k_k_method (choises: 1, 2, 3) you can chose the method for calculating the matrix elements. 1 will first calculate the matrix elements (k, k’)!!! for the atomic basis before transforming to the KS-basis (biggest array size: n_basis*n_k_points*n_k_points). 2 will calculate the matrix elements (k) for the atomic basis and transform to the KS-basis for all (k’) (n_k_points times) (biggest array size: n_basis*n_k_points). 3 will calculate the matrix elements in the atomic basis and transform to the KS-basis for all (k, k’) (n_k_points*n_k_points times) (biggest array size: n_basis). 1 needs the most memory, 3 does the most (repetitive) calculations.
The dipole matrix elements for all (k,k’)-points are written into one container file (dipmat_k.h5). To use this functionality, FHI-aims has to be compiled with the external hdf5 module (see Sec. H.2 for details).

Tag: compute_kubo_greenwood(control.in)

Usage: compute_kubo_greenwood width FD_Temp $E_{min}$ $E_{max}$ $\omega_{min}$ $\omega_{max}$ $n_{\omega}$ i j
Purpose: Calculate and output the i j component (choices: (x, y or z) or (a a)) of the the Kubo-Greenwood transport properties optical conductivity $\sigma(\omega)$ and Seebeck coefficient S in units $(\Omega\cdot cm)^{-1}$ and $\mu V/K$ , respectively. The (a a) choice triggers the calculation of the average of the three diagonal elements $L_{ij,xx}$ , $\sigma_{ij,yy}$ and $\sigma_{ij,zz}$ . As before, [ $E_{min}$ , $E_{max}$ ] (in eV) determine the energy window (relative to the internal zero) within which the momentum-matrix elements are considered. The Delta-distribution is represented by a Gaussian function with width width (in eV) and the Fermi occupations are calculated at an electronic temperature FD_Temp (in eV). The spectra contain $n_{\omega}$ $\omega$ -values in the interval [ $\omega_{min}$ , $\omega_{max}$ ] (in eV). The output file is named dielectric_i_j<Fermi-level>.out and consists of the dielectric function, the optical conductivity and the Seebeck coefficient. Usually it is enough to include only a few states below and above the fermi level in the energy window [ $E_{min}$ , $E_{max}$ ].
The calculation is very sensitive to the k-point grid, an extremely high number of k-points might be needed for convergence, especially for metals.

Tag: greenwood_method(control.in)

Usage: greenwood_method method
Purpose: Switches between the use of sparse or full transition matrices. (Choices: sparse or full). As of now, the keyword has to be set in case the compute_kubo_greenwood is used. The use of sparse matrices is recommended. (Keyword will be deprecated in the near future)[1.0ex]

Tag: calculate_plasma_frequency(control.in)

Usage: calculate_plasma_frequency
Purpose: Activates the calculation of the spin-resolved plasma frequency tensor. Calculations are performed as a post-SCF step. This keyword must be used in combination with dos_kgrid_factors.

When using this module in your work, please cite the following reference:

•

M. K. Pogodaeva, S. V. Levchenko, and V. P. Drachev, “Spin-dependent plasma frequency from all-electron ab initio calculations including spin-orbit coupling,” Phys. Rev. B 107, 045113 (2023), https://doi.org/10.1103/PhysRevB.107.045113.

Usage example

To calculate the plasma frequency, add the following lines to your control.in:

    dos_kgrid_factors n_x n_y n_z
    calculate_plasma_frequency

Here, n_x, n_y, n_z are factors specifying how much denser the new k-grid is compared to the original SCF cycle. If no denser grid is required, use:

    dos_kgrid_factors 1 1 1

Optionally, you can specify to output the density of states computed on the k-grid defined by dos_kgrid_factors at no additional computational expense, via the keyword output dos. Please see the sections of the manual where these are defined for further details.

If you have to encounter SOC, simply add include_spin_orbit keyword. The code will produce both spin-perturbed and spin-unperturbed plasma frequency.

If you need spin-resolved plasma frequency, specify spincollinear and the plasma frequency tensor for both spin channels will be printed in the standard output.

Note, that you do not need compute_dielectric keyword to run plasma frequency calculations.