Spectroscopic Properties

Time-dependent DFT

In this section, the time-dependent density functional theory implementation in BAND is described. How to do a response calculation in BAND (which input keys to use), can be found in the list of keywords, see the key Response. The TDDFT module enables the calculation of real and imaginary parts of the material property tensor \(\chi_e(\omega)\) called the electric susceptibility, and the macroscopic dielectric function \(\epsilon_e(\omega)\) These are mutually related,

\[\epsilon_e(\omega) = 1 + 4 \pi \chi_e(\omega)\]

In general \(\chi_e(\omega)\) and \(\epsilon_e(\omega)\) are tensors, which, however, simplify to scalars in isotropic systems. The above formula is valid in the case of insulators and semiconductors, where the bands are either fully occupied or fully unoccupied. It is defined as the interband part of the dielectric function and it is due to transitions from occupied bands to unoccupied bands. Also in the metallic case similar expressions can be found, the main difference is that there now is also an intraband contribution. Some examples are available in the $ADFHOME/examples/band directory and are discussed in the Examples document.

References

The three related Ph.D. theses, due to F. Kootstra (on TD-DFT for insulators), P. Romaniello (on TD-CDFT for metals), and A. Berger (on the Vignale-Kohn functional in extended systems) contain much background information, and can be downloaded from the SCM website.

The most relevant publications on this topic due to the former “Groningen” group of P.L. de Boeij are [22], [23], [24], [25].

Response key

Response

Perform a time-dependent DFT calculation to obtain real and imaginary parts of frequency-dependent dielectric function.

Response
   {nfreq   ival1}
   {strtfr  rval1}
   {endfr   rval2}
   {cnvi    rval3}
   {cnvj    rval4}
   {ebndtl  rval5}
   {shift   rval6}
   {isz     ival2}
   {iyxc    ival3}
   {static}
   {newvk}
   {cnt}
   {qv}
   {Berger2015}
End

Omitting the specific options in the Response block will cause default setting to be used during the calculation, as given above.

nfreq

(Default: 5) the number of frequencies in a.u. for which a TDDFT calculation is performed when calculating the dielectric function \(\epsilon_e(\omega)\) of a system.

strtfr

(Default: 0d0) is the start frequency in a.u. of the frequency range over which the dielectric function is calculated.

endfr

(Default: 1d-2) is the end frequency in a.u. of the frequency range over which the dielectric function is calculated.

cnvi

(Default: 1d-3) the first convergence criterion for the change in the fitcoefficients for the fitfunctions, when fitting the density.

cnvj

(Default: 1d-3) the second convergence criterion for the change in the fitcoefficients for the fitfunctions, when fitting the density.

ebndtl

(Default: 1d-3) the energy band tolerance, for determination which routines to use for calculating the numerical integration weights, when the energy band posses no or to less dispersion.

shift

(Default: 0d0) shift (in a.u.) of the virtual crystal orbitals.

isz

(Default: 0) integer indicating whether or not scalar zeroth order relativistic effects are included in the TDDFT calculation. 0 = relativistic effects are not included, 1 = relativistic effects are included.

iyxc

(Default: 0) integer for printing yxc-tensor (JCP 115, 1995 (2001)). 0 = not printed, 1 = printed.

static

An alternative method that allows an analytic evaluation of the static response (Normally the static component is approximated by a finite small value). This option should only be used for non-relativistic calculations on insulators, and it has no effect on metals. Experience shows that KSPACE convergence can be slower. That is why it is not the default.

newvk

Use the slightly modified version of the VK kernel.[58] When using this option one uses effectively the static option, even for metals, so one should check carefully the convergence with the KSPACE parameter. An example response block for a vignale kohn calculation looks like

qv/cnt

Use the QV or CNT parametrization for the longitudinal and transverse kernels of the xc-kernel of the homogeneous electron gas. Use this in conjunction with the newvk option.[60,61]

Response
   {...}
   newvk
   iyxc    1
   cnt
END

Berger2015

Use the parameter-free polarization functional by A. Berger.[59] This is possible for 3D insulators and metals, but not in combination with relativistic approximations.

Response
   {...}
   Berger2015
END

Limitations

The method has not been implemented for slabs, so it can only be used for 1D and 3D systems.

References

[2] F. Kootstra, P. L. de Boeij, and J. G. Snijders, Phys. Rev. B 62, 7071 (2000).

[3] F. Kootstra, P. L. de Boeij, H. Aissa, and J. G. Snijders, J. Chem. Phys. 114, 1860 (2001).

[4] P. L. de Boeij, F. Kootstra, and J. G. Snijders, Int. J. Quantum Chem. 85, 449 (2001).

[5] P. L. de Boeij, F. Kootstra, J. A. Berger, R. van Leeuwen, and J. G. Snijders, J. Chem. Phys. 115, 1995 (2001).

[6] F. Kootstra, P. L. de Boeij, R. van Leeuwen, and J. G. Snijders, Festschrift in honour of R. G. Parr, Editor K. D. Sen, accepted.

[7] F. Kootstra, P. L. de Boeij, and J. G. Snijders, J. Chem. Phys. to be submitted.

[8] F. Kootstra, P. L. de Boeij, R. van Leeuwen, and J. G. Snijders, to be submitted.

[9] F. Kootstra, Ph.D. thesis, Rijksuniversiteit Groningen, Groningen (2001).

Time-dependent DFT for metals

For metals the interband part of the dielectric function is due to transitions from (partially)-occupied bands to (partially)-unoccupied bands. Now there is also a term, which is called the intraband part of the dielectric function, which is due to transitions within the same partially-occupied band. The macroscopic dielectric function \(\epsilon_e(\omega)\) is now calculated as

\[\epsilon_e(\omega) = 1 + 4 \pi \chi_e(\omega) - 4 \pi i \sigma_e(\omega) / \omega\]

Convergence and reproducibility: For TD-DFT calculations on metals a dense sampling of reciprocal space is required, i.e. the results converge slowly with the KSPACE parameter. The dielectric function might even not reproduce well across different machines. (Usually results that are not yet converged with the KSPACE parameter, like the energy, are reproducible across different platforms). Nevertheless, when the KSPACE parameter is chosen sufficiently high, the same result will be obtained on all machines. For instance for Cu the machine dependence was less than 0.01 for the dielectric function with KSPACE=11. In short: check for the convergence of the dielectric function with respect to the KSPACE parameter.

Frequency dependent kernel

It is known that the exact Vignale Kohn kernel greatly improves the static polarizabilities of infinite polymers and nanotubes (see JCP 123 174910), but gives bad results for the optical spectra of semiconductors and metals. For the low frequency part one needs a frequency dependent kernel, because Drude-like tails are completely absent in the ALDA. With a modified Vignale Kohn kernel, neglecting \(\mu_{xc}\) so that it reduces to the ALDA form in the static limit (see PRB 74 245117) much better results can be obtained. BAND currently only supports the modified VK kernel in either the QV or CNT parametrization, and it should only be used for metals.

EELS

Once the macroscopic dielectric function is known it is possible to calculate the electron energy loss function (EELS). In transmission electron energy loss spectroscopy one studies the inelastic scattering of a beam of high energy electrons by a target. The scattering rates obtained in these experiments is related to the dynamical structure factor \(S(q,\omega)\) [1]. In the special case with wavevector \(q=0\), \(S(q,\omega)\) is related to the longitudinal macroscopic dielectric function. This is the long-wave limit of EELS. For isotropic system the dielectric function is simply a scalar ( \(1/3 \text{Tr} (\epsilon_e(\omega))\) ). In this case the long-wave limit of the electron energy loss function assumes the trivial form

\[\lim_{q \rightarrow 0} 2 \pi \frac{S(q,\omega)}{q^2 V} = \frac{\epsilon_2}{\epsilon_1^2 + \epsilon_2^2}\]

with \(\epsilon_1\) and \(\epsilon_2\), respectively, the real and imaginary part of the dielectric function.

[1] S. E. Schnatterly, in Solid State Physics Vol.34, edited by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic Press, Inc., New York, 1979). [2] P. Romaniello, and P. L. de Boeij, Phys. Rev. B (accepted).

ESR

BAND is able to calculate electron paramagnetic resonance parameters of paramagnetic defects in solids: hyperfine A-tensor and the Zeeman g-tensor.

The implementation of EPR parameters in BAND are described in the publications by Kadantsev and coworkers [20] and [21].

Hyperfine A-tensor

The A-tensor is implemented within the non-relativistic and scalar relativistic spin-polarized Kohn-Sham scheme. The A-tensor calculation is invoked by block

ATENSOR
END

Note: the Unrestricted keyword should be present.

Two methods are used for A-tensor calculation.

Method 1 involves gradient of spin-polarization density and integration by parts. The isotropic component of A-tensor is obtained through integration, in a “nonlocal fashion”.

In Method 2, the A-tensor is computed from spin-polarization density. Method 2 does not relies on the integration by parts. The isotropic component is obtained in a “local fashion” from the value of spin-polarization density on the grid points near the nuclei.

The user should be aware that numerical integration in A- and g-tensor routines is carried out over Wigner-Seitz (WS) cell, and, therefore, to obtain a meaningful result, the defect in question should lie at or, very close to, WS cell origin. This might require, on the user’s part, some modification of the input geometry.

It also might happen that the size of the WS cell is not large enough for the adequate description of the paramagnetic defect in question. In this case, Method 1, which relies on the integration by parts and assumes that the spin-polarization density is localized inside the WS cell will fail. For the same reason, We recommend that the user removes diffuse basis set functions that describe the defect subsystem.

Finally, we note that the final result for A-tensor as presented by BAND is not scaled by the nuclear spin (as it is done in ADF) and the user is responsible for making necessary adjustments.

g-tensor

The calculation of Zeeman g-tensor is invoked with block

ESR
END

The Zeeman g-tensor is implemented using two-component approach of Van Lenthe and co-workers in which the g-tensor is computed from a pair of spinors related to each other by time-reversal symmetry.

Note: The keyword Relativistic zora spin should be present to invoke calculations with spin-orbital coupling. The user also has to specify

Kspace 1

(\(\Gamma\)-point calculation). The g-tensor is then computed from the HOMO spinor at the \(\Gamma\) point. In the output, the user can find two-contributions to g-tensor: one that stems from \(K_\sigma\) operator and a second one, that stems from orbital angular momentum. By default, GIAO and spin-Zeeman corrections are not included, from our experience these corrections are quite small.

Electric Field Gradient (EFG)

EFG (block type)

The electronic charge density causes an electric field, and the gradient of this field couples with the nuclear quadrupole moment, that some (non-spherical) nuclei have and can be measured by several spectroscopic techniques. The EFG tensor is the second derivative of the Coulomb potential at the nuclei. For each atom it is a 3x3 symmetric and traceless matrix. Diagonalization of this matrix gives three eigenvalues, which are usually ordered by their decreasing absolute size and denoted as \(|V_{zz}|\) , \(|V_{yy}|\), \(|V_{xx}|\). The result is summarized by the largest eigenvalue and the asymmetry parameter.

This options honors the SelectedAtoms key, in which case only the EFG will be calculated for the selected atoms.

NMR

With the NMR option the shielding tensor is calculated. There are essentially two methods: the super cell method and the single-dipole method.

1. The super cell method is according to the implementation by Skachkov et al..[37] The symmetry will automatically be set to NOSYM. The unit cell should not be chosen too small.
2. The other method is the single-dipole method. In principle one can now use the primitive cell.[48] In practice also this method needs to be converged with super cell size. However, depending on the system the required super cell may be much smaller. At a given super cell size this method is more expensive than the super cell method.

NMR (block type)
   {SuperCell [true | false]}
End

SuperCell

(Default: true) This is the switch between the two methods, either the super cell (true), or the single-dipole method (false).

This options honors the SelectedAtoms key, in which case only the NMR properties will be calculated for the selected atoms only.