DFT: LCAO

Introduction

QuantumATK can model the electronic properties of closed and open quantum systems within the framework of density functional theory (DFT) using numerical LCAO basis sets (linear combination of atomic orbitals). For periodic systems, QuantumATK can also use the plane-wave basis set for DFT calculations, as discussed in DFT: Plane Wave.

The key parameter in the self-consistent calculation of the Kohn–Sham equations is the density matrix, which defines the electron density. For open systems, the density matrix is calculated using non-equilibrium Green’s functions (NEGFs), see NEGF: Device Calculators, while for closed or periodic systems it is calculated by diagonalization of the Kohn–Sham Hamiltonian, as described in the present chapter. The electron density then sets up an effective potential, which is given by the Hartree, exchange-correlation, and external potentials. Knowing the effective potential allows us to obtain the Kohn–Sham Hamiltonian.

The next section describes the mathematical formalism behind the DFT-LCAO model.

Background information

The LCAOCalculator provides a description of electronic structure using DFT and norm-conserving pseudopotentials. The method is based on an expansion of the single-particle wave functions in a basis of numerical atomic orbitals with compact support. In this section, we describe the mathematical formalism of the DFT-LCAO calculator, which closely follows the work of Soler et al. [1].

Kohn–Sham Hamiltonian

Within density functional theory, the many-body electronic structure of the system is described in terms of the one-electron Kohn–Sham Hamiltonian:

\[\hat{H}_\mathrm{1el} =-\frac{\hbar^{2}}{2m} \nabla^2 + V^{\mathrm{eff}}[n](\mathbf{r}).\]

In this equation, the first term is the kinetic energy of the electron, while the second term (the effective potential) is the potential energy of the electron moving in the mean field created by the other electrons as well as in any external potential field, e.g. the electrostatic potential of ions or any other external field. The electrons are described in terms of the total electron density, \(n=n(\mathbf{r})\).

The electron density is discussed in detail in the section Electron density, and the effective potential in the section Effective potential.

Solving the Kohn–Sham equations by means of a basis set expansion

We calculate the one-electron eigenfunctions of the Kohn–Sham Hamiltonian, \(\psi_{\alpha}\), by solving the one-electron Schrödinger equation,

\[\hat{H}_{\mathrm{1el}} \psi_{\alpha}({\mathbf{r}}) = \varepsilon_{\alpha} \psi_{\alpha}({\mathbf{r}}).\]

This differential equation is called the Kohn–Sham equation within DFT. To solve it, we expand the eigenfunctions \(\psi_{\alpha}({\mathbf{r}})\) in a set of basis functions, \(\phi_i\):

\[\psi_{\alpha}(\mathbf{r}) = \sum_i c_{\alpha i} \phi_i(\mathbf{r}).\]

This allows us to represent the differential equation as a matrix equation for determining the expansion coefficients, \(c_{\alpha i}\):

\[\sum_{j } H_{ij} c_{\alpha j} = \varepsilon_\alpha \sum_{j } S_{ij} c_{\alpha j},\]

where the Hamiltonian matrix, \(H_{ij} = \langle \phi_i | \hat{H}_{\mathrm{1el}} | \phi_j \rangle\), and overlap matrix \(S_{ij} = \langle \phi_i | \phi_j \rangle\) are given by the multiple integrals with respect to the electron coordinates.

Electron density

The electron density of the many-electron system is given by the occupied eigenstates of the Kohn–Sham Hamiltonian:

\[n(\mathbf{r}) = \sum_{\alpha} f_\alpha |\psi_\alpha(\mathbf{r})|^2,\]

where \(f_\alpha\) is the occupation of the level denoted by \(\alpha\). For finite temperature calculation the occupations are determined by the Fermi-Dirac distribution \(f_\alpha = \frac{1}{1 + e^{(\epsilon_\alpha - \epsilon_F)/kT}}\) but other smooth distributions may be introduced in order to help speed convergence (see Occupation Methods).

The electron density can then be expressed in terms of the density matrix:

\[n(\mathbf{r}) = \sum_{ij} D_{ij} \phi_i(\mathbf{r}) \phi_j(\mathbf{r}),\]

where the density matrix is given by the basis set expansion coefficients \(c_{\alpha i}\):

\[D_{ij} = \sum_{\alpha} f_\alpha c_{\alpha i}^* c_{\alpha j}.\]

Note

For open systems (DeviceConfiguration), the density matrix is calculated using non-equilibrium Green’s functions, see NEGF: Device Calculators.

Electron difference density

It is often convenient to compare the electron density of the many-body system to a superposition of individual atom-based electron densities, \(n^{\mathrm{atom}}(\mathbf{r}-\mathbf{R}_\mu)\), where \(\mathbf{R}_\mu\) is the position of atom \(\mu\) in the many-body system:

\[\Delta n(\mathbf{r}) = n(\mathbf{r}) - \sum_{\mu} n^{\mathrm{atom}}(\mathbf{r} - \mathbf{R}_\mu).\]

\(\Delta n\) is called the electron difference density, and it is calculated using the ElectronDifferenceDensity analysis object.

Effective potential

The effective potential, \(V^\mathrm{eff}[n]\), has three contributions:

\[V^{\mathrm{eff}}[n] = V^{H}[n] + V^\mathrm{xc}[n] + V^\mathrm{ext}.\]

The first two terms are due to electron–electron interactions, which depend on the electron density. The first term, \(V^{H}[n]\), is the Hartree potential due to the mean-field electrostatic interaction between the electrons, while the second term, \(V^\mathrm{xc}[n]\), is the exchange-correlation potential, which arises from the quantum mechanical nature of the electrons.

The potential \(V^\mathrm{ext}\) represents any other electrostatic fields in the system. It can be separated into two contributions; the electrostatic potential of ions (given by norm-conserving pseudopotentials) and external electrostatic fields (given by one or more external sources).

External potential

The external potential is given by the pseudopotentials and an external electrostatic field. Such a field may arise from the inclusion of metallic gates:

\[V^\mathrm{ext} = \sum_\mu V^\mathrm{pseudo}_\mu + V^\mathrm{gate}.\]

The pseudopotential has two contributions:

\[V^\mathrm{pseudo}=V^\mathrm{local}+\sum_{n,n'}|\chi_n \rangle B_{n,n'}\langle \chi_{n'} |,\]

where the first and second terms correspond to the local and nonlocal part of the pseudopotential, respectively.

The term \(V^\mathrm{gate}\) is the electrostatic potential due to external gates, calculated with zero electron density. This term will be returned by the Analysis object ExternalPotential and is calculated using the MultigridSolver, DirectSolver and ParallelConjugateGradientSolver Poisson solvers. The Analysis object ElectrostaticDifferencePotential also includes this term, among others.

Hartree potential

Using the electron density, we can calculate the classical electrostatic potential, the so-called Hartree potential. The actual calculation of the Hartree potential is described in detail in the section Poisson solvers.

Exchange-correlation potential

In the DFT method, the quantum mechanical part of the electron–electron interaction is approximated by the exchange-correlation term, and a large number of different approximate exchange-correlation density functionals exist. ATK supports many of these, see the section Exchange-correlation energy.

The exchange-correlation potential is defined as the functional derivative of the exchange-correlation energy with respect to the electron density, which corresponds to a mean-field quantum mechanical interaction potential between the electrons:

\[V^\mathrm{xc}[n](\mathbf{r}) = \frac{\delta E^\mathrm{XC}} {\delta n}(\mathbf{r}).\]

Total energy and forces

The DFT total energy of a many-electron system is a functional of the electron density, \(n\):

\[E[n] = T[n] + E^\mathrm{xc}[n] + E^{H}[n] + E^\mathrm{ext}[n],\]

where \(T[n]\) is the kinetic energy of a non-interacting electron gas with density \(n\), \(E^\mathrm{xc}[n]\) the exchange-correlation energy, \(E^{H}[n]\) the Hartree potential energy, and \(E^\mathrm{ext}[n]\) the interaction energy of the electrons in the electrostatic field created by ions and other external sources.

The electron kinetic energy may be defined as

\[T[n] = \sum_{\alpha} f_\alpha \langle \psi_{\alpha}| \frac{-\hbar^{2}}{2m} \nabla^{2} | \psi_{\alpha} \rangle.\]

The total energy is calculated using TotalEnergy.

First-principles forces are calculated by differentiating the total energy with respect to the ionic coordinates of atom \(i\) at position \(\mathbf{R}_i\):

\[\mathbf{F}_i = -\frac{d E[n]}{d \mathbf{R}_i}.\]

Pseudopotentials

QuantumATK uses norm-conserving pseudopotentials and PAW potentials, and is shipped with a database for the entire periodic table. See the full list here: Pseudopotentials. This database is reviewed and updated on a regular basis to provide increasingly accurate and general-purpose pseudopotentials and PAW potentials.

QuantumATK uses the unified pseudopotential format (upf) defined by the Quantum ESPRESSO consortium. From their website, one may download tools to convert several different pseudopotential formats into a upf file.

LCAO basis set

The eigenfunctions of the Kohn–Sham Hamiltonian can be expanded in a Linear Combination of Atomic Orbitals (LCAO’s):

\[\phi_{nlm}(\mathbf{r}) = R_{nl}(r) Y_{lm}(\hat{\mathbf{r}}),\]

where \(Y_{lm}\) are spherical harmonics, and \(R_{nl}\) are radial functions with compact support, being exactly zero outside a confinement radius.

Note

The basis set functions have a finite range, but the interaction range of the Hamiltonian is larger than that of the basis set. Theoretically, the Hamiltonian interaction range may exceed twice the basis set range, but is usually smaller in practice.

The basis orbitals have a number of parameters that determine the shape of the orbitals. It is possible to assemble the basis orbitals into your own basis set through the use of the BasisSet keyword.

QuantumATK comes with a number of pre-built basis sets for each chemical element and for each type of pseudopotential. See the full list here: LCAO basis sets

Exchange-correlation energy

The QuantumATK package use the libxc library for providing a large suite of exchange-correlation functionals. For each functional it is possible to add a Hubbard correction, as described in the section XC+U mean-field Hubbard term.

The main families of exchange-correlation functionals supported in QuantumATK are the Local-density approximation (LDA), the Generalized-gradient approximation (GGA), and the Meta-GGA. Each functional comes in four spin variants:

  • non-polarized,

  • spin-polarized (collinear),

  • noncollinear,

  • noncollinear with spin-orbit coupling.

Important

The choice of spin variant for the exchange-correlation functional determines the spin type used in all of the DFT-LCAO calculation.

Furthermore, a Hubbard correction can be added to all exchange-correlation variants, regardless of spin type, see XC+U mean-field Hubbard term. Several standard GGA functionals can also be extended with the DFT-D2 and DFT-D3 dispersion corrections by Grimme and co-workers.

Initial spin

Both collinear and noncollinear calculations may have local minima of the total energy that correspond to different spin configurations. It is therefore important to prepare the system in the correct initial spin configuration, such that the DFT self-consistent calculation will end up in the lowest energy spin configuration. This is done through the use of the methods InitialSpin and RandomSpin. The initial spin direction on each atom is specified in physical spherical coordinates \((r,\theta,\phi)\), where \(\theta\) is the angle with the z-axis, and \(\phi\) the polar angle in the x-y plane relative to the x-axis. The collinear case is \(\theta=0\) Radians or \(\theta=\pi\) Radians.

Local-density approximation (LDA)

In the LDA approximation, the exchange-correlation energy is taken as a functional of the local electron density,

\[E^\mathrm{LDA}[n] = \int n(\mathbf{r}) \varepsilon^\mathrm{LDA}(n(\mathbf{r})) d\mathbf{r},\]

where \(\varepsilon^\mathrm{LDA}(n(\mathbf{r}))\) is the exchange-correlation energy density of a homogeneous electron gas with density \(n(\mathbf{r})\).

It is possible to derive an exact, analytical expression for the exchange energy of the homogeneous electron gas, the so-called Dirac–Bloch exchange energy, which is used for all the LDA functionals. The correlation energy of the electron gas cannot be calculated exactly, so a number of different approximations have been proposed over the years. Most of them give almost identical results; the most commonly used variant is the parametrization by Perdew and Zunger, which is also the default LDA functional in DFT: LCAO.

For a list of the parametrizations for the LDA correlation functional implemented in QuantumATK, see the documentation of the ExchangeCorrelation class. In particular, the section Abbreviations explains how to select a particular exchange-correlation functional. See also Table 18, Table 21, and Table 22.

Generalized-gradient approximation (GGA)

The GGA functionals is a large family of semi-local approximations for the exchange-correlation energy, where the functional depends on both the local value and the local gradient of the electron density,

\[E^\mathrm{GGA}[n] = \int n(\mathbf{r}) \varepsilon^\mathrm{GGA}(n(\mathbf{r}), \nabla n(\mathbf{r})) d\mathbf{r}.\]

For a list of the GGA exchange-correlation functionals available in QuantumATK, see the documentation of the ExchangeCorrelation class. In particular, the section Abbreviations explains how to select a particular exchange-correlation functional. See also Table 19, Table 21, and Table 22.

Meta-GGA

The meta-GGA (MGGA) family of functionals extend the GGA approximation by additionally depending on one or both of the following quantities: the Laplacian of the density

\[\nabla^2 n(\mathbf{r})\]

and the so-called kinetic energy density

\[\frac{1}{2} \sum_\alpha f_\alpha \left | \nabla \psi_\alpha (\mathbf{r}) \right |^2.\]

For a list of the MGGA exchange-correlation functionals available in QuantumATK, see the documentation of the ExchangeCorrelation class. In particular, the section Abbreviations explains how to select a particular exchange-correlation functional. See also Table 20, Table 21, and Table 22.

Different formulations of the MGGA can produce functionals with different uses. The SCAN functional has been found to give improved energetics over the LDA and GGA in several cases; however, band gaps are similar.

Conversely, the TB09 functional can provide a much more accurate description of band gaps than ordinary LDA and GGA. In fact, the accuracy of semiconductor band gaps obtained with TB09 are often comparable to those from calculations using GW or hybrid functionals, which have a significantly higher computational cost. It is even possible to tune the TB09 method such that it works optimally on both sides of an interface (see ExchangeCorrelation). Alternatively, one can use the Hubbard U model within the LDA or GGA with the U parameters for the semiconductor elements fitted against TB09 results.

However, TB09 is not meant for energetics, and, hence, geometry optimizations are disabled with this functional. For this, the GGA or other MGGA functionals should be used.

XC+U mean-field Hubbard term

The local approximations to the exchange-correlation energy have a number of shortcomings. Two of these are of particular interest:

  • Self-interaction: The electron is formally allowed to interact with itself. This can prevent electrons from localizing properly.

  • Excited states: The LDA and GGA description of conduction-band energy levels is often poor, so band gaps are often too low.

The mean field Hubbard correction by Dudarev et al. [2] and Cococcioni et al. [3], often denoted XC+U, DFT+U, LDA+U, or GGA+U, is a semi-empirical correction which attempts to improve on these deficiencies of the local exchange-correlation functionals by adding an extra term to the exchange-correlation functional:

\[E_{U} = \frac{1}{2} \sum_\mu U_\mu (n_\mu - n_\mu^2) .\]

In this equation, \(n_\mu\) is the projection onto an atomic shell and \(U_\mu\) is the Hubbard U for that shell. The \(E_{U}\) energy term is zero for a fully occupied or unoccupied shell, while positive for a fractionally occupied shell.

The energy is thereby lowered if states become fully occupied or empty. This may happen if the energy levels move away from the Fermi Level, increasing the band gap, or if the broadening of the states is decreased, the electrons are then more localized. Thus, the Hubbard U method improves on the deficiencies of the local exchange-correlation functionals listed above.

The value of \(U\) is often used as an empirical parameter, which is varied in order to improve the comparison between DFT and experimental data. However, it has also been suggested that the value of U may be obtained by minimizing the total energy of the system [3].

Some caution must be taken when using the XC+U correction, since for large U values the electron density may have several local minima where some of them are unphysical, and it is often necessary to select an anisotropic initial electron state to reach the correct local minimum.

XC+U implementations in QuantumATK

The XC+U implementation in QuantumATK comes in four main variants; a standard local-orbital formulation, so-called Onsite representation, and the Dual representation introduced by Han et al. [4]. The difference between the two implementations is in the definition of the local occupation matrix. The occupations may be summed over a single orbital or over an entire angular momentum shell, the latter corresponds to the OnsiteShell and DualShell representations.

Onsite representation

In the onsite representation, the local occupation matrix is given by

\[n_{\mu, m m'}^\sigma = D_{\mu m, \mu m' }^\sigma,\]

where \(D\) is the density matrix, \(\mu\) a basis orbital index, and \(\sigma\) a spin index.

OnsiteShell representation

In the onsite shell representation, the occupation is obtained by summing together all the basis functions in each angular momentum shell, \(l\):

\[n_{l, m m'}^\sigma = \sum_{\mu \in l} n_{\mu, m m'}^\sigma.\]
Dual representation

In the dual representation, the local occupation matrix is given by

\[n_{\mu, m m'}^\sigma = \sum_{i, n} \left[ S_{\mu m, i n} D_{i n, \mu m'}^\sigma + D_{\mu m, i n }^\sigma S_{ i n, \mu m'} \right],\]

where \(S\) is the overlap matrix.

DualShell representation

Similar to the onsite shell representation, the dual shell occupation is obtained by summing up all dual occupations in each angular momentum shell, \(l\):

\[n_{l, m m'}^\sigma = \sum_{\mu \in l} n_{\mu, m m'}^\sigma.\]

Which model to use?

The onsite representation corresponds to projecting into a local orbital, and since the local orbitals are non-orthogonal, there are no sum rules for the occupation matrix and it is not guaranteed that a fully occupied shell has occupation 1.

In the dual representation, the occupation matrix has the form of a Mulliken population, and this form has the advantage that sum rules apply, i.e. the trace of the occupation matrix sums up to the total number of electrons.

The dual representation corresponds to projecting into orthogonalized orbitals. Such orbitals may have non-zero weights also on neighboring centers, and our experience shows that this kind of non-locality can result in non-physical results in some cases. The default in QuantumATK is therefore to use the onsite representation.

The onsite shell and dual shell representations often give a better description with multiple-zeta basis sets compared to the onsite and dual representation, where the occupation is obtained by projecting into single basis functions.

Input format for XC+U

The Hubbard U approximation in the onsite representation may be selected through keywords of the types LDAU.XCTYPE, LSDAU.XCTYPE, GGAU.XCTYP, and SGGAU.XCTYPE.

The keyword LSDAU.PZ gives the default LSDA functional with the default Hubbard U method. The specification

exchange_correlation = LSDAU.PZ

is therefore identical to

exchange_correlation = ExchangeCorrelation(
    exchange=DiracBloch,
    correlation=PerdewZunger,
    hubbard_term=Onsite,
    number_of_spins=2,
    )

Use this form with the option hubbard_term=Dual to select the dual representation.

The value of the U term is specified through the basis set:

basis_set = [LDABasis.Nickel_SingleZeta(hubbard_u=[4.6, 0.0]*eV,
                                        filling_method=Anisotropic)]

which will give a single-zeta basis set for nickel, consisting of a 3d and a 4s orbital, where the 3d orbital has U=4.6 eV. The 3d orbitals are filled anisotropically, i.e, there is 1 electron in in orbitals with m=-2,-1,0,1, leaving m=2 empty. The alternative is filling_method=SphericalSymmetric, which puts 4/5 electron in each orbital.

DFT-1/2 method

The DFT-1/2 method (often also denoted LDA-1/2 or GGA-1/2) is a semi-empirical approach to correct the self-interaction error in local and semi-local exchange-correlation functionals for extended systems. As such, it can be viewed as an alternative to the XC+U mean-field Hubbard term with broadly the same aims: to improve the description of conduction-band energy levels and band gaps.

This method, introduced by Ferreira et al. [5], is based on the much older Slater half-occupation scheme for molecules [6] (also known as the transition state method). Slater’s original method consists in carrying out a self-consistent calculation with half an electron removed from the system, and taking the eigenvalue of the half-filled state as an estimate for the ionization energy. This approach was later formalized in the theorem by Janak [7]:

\[\frac{\partial E}{\partial f_\alpha} = \varepsilon_\alpha \left ( f_\alpha \right ),\]

where \(E\) is the total energy of the system, \(f_\alpha\) is the occupation of state \(\alpha\) (between 0 and 1), and \(\varepsilon_\alpha\) is the eigenvalue of the state. The success of the half-occupation scheme is based on the fact that \(\varepsilon_\alpha \left ( f_\alpha \right )\) is known to be almost precisely linear for many cases, which makes the relationship with the ionization energy exact.

The DFT-1/2 method makes use of the theoretical insights from the half-occupation scheme and to tackle the fundamental problem of the self-interaction error for the case of extended systems with a band gap. The method attempts to correct for this by defining an atomic self-energy potential which cancels the electron-hole self-interaction energy. This potential is calculated for atomic sites in the system, and is defined as the difference between the potential of the neutral atom and that of a charged ion resulting from the removal of a fraction of its charge, between 0 and 1 electrons. The total self-energy potential is the sum of these atomic potentials.

The addition of the DFT-1/2 self-energy potential to the DFT Hamiltonian has been found to greatly improve band gaps for a wide range of semiconducting and insulating systems [8].

It is important to note that the method is not entirely free of empirical parameters. Much like the value of \(U\) for the Hubbard U method, there are two values which much be fixed for every species in the system: the fractional charge \(f_\alpha\) removed from the neutral atom, and a cutoff radius \(r^\mathrm{cut}\) beyond which to trim the atomic self-energy potential. The latter is needed to avoid excessive overlap of the potentials from different sites in the crystal.

The original authors of the method suggest to fix \(r^\mathrm{cut}\) variationally, by choosing the value which maximizes the band gap [5] [8].

The choice of \(f_\alpha\) is less well-defined: although a value of 0.5 is standard (hence the name of the method), this is known not to be the optimal choice for various materials (e.g., silicon). It may therefore be appropriate to treat it as an empirical parameter, to be varied by comparison with experiment.

It is also important to note that not all species in the system necessarily require the DFT-1/2 correction; it is generally advisable only to add this to the anionic species, and leave the cationic species as normal [5] [8].

Shell-DFT-1/2

Following the paper by Xue et al. [9] QuantumATK also supports calculations with the slightly modified shell-DFT-1/2 method. The difference between normal DFT-1/2 and shell-DFT-1/2 is that in the latter, the self-energy potential is added in a shell around the atom, and not in a sphere as in normal DFT-1/2. The shell is specified by an additional inner radius parameter in addition to the cutoff radius \(r^\mathrm{cut}\).

Input format for DFT-1/2

The DFT-1/2 method may be selected by specifying one of the pre-defined exchange-correlation functionals: LDAHalf, LSDAHalf, NCLDAHalf, SOLDAHalf, GGAHalf, SGGAHalf, NCGGAHalf, SOGGAHalf.

The keyword LDAHalf.PZ gives the default LDA functional with the DFT-1/2 correction. The specification

exchange_correlation = LDAHalf.PZ

is therefore identical to

exchange_correlation = ExchangeCorrelation(
    exchange=DiracBloch,
    correlation=PerdewZunger,
    dft_half_enabled=True
    )

In general, the keyword dft_half_enabled=True can be specified for any user-defined exchange-correlation functional.

The value of the DFT-1/2 parameters are specified through the basis set:

dft_half_parameters = DFTHalfParameters(
    element=Arsenic,
    fractional_charge=[0.3, 0.0],
    cutoff_radius=4.0*Bohr,
    inner_radius=0.0*Bohr)

basis_set = [
    LDABasis.Arsenic_DoubleZetaPolarized(
        dft_half_parameters=dft_half_parameters),
    LDABasis.Gallium_DoubleZetaPolarized(
        dft_half_parameters=Disabled)
]

The above example shows a typical case for gallium arsenide, in which the anion (As) is given a DFT-1/2 self-energy potential calculated with \(f_\alpha = 0.3\) and \(r^\mathrm{cut} = 4~a_0\), while the cation (Ga) has its DFT-1/2 correction turned off with the keyword Disabled.

If dft_half_parameters is not specified for a species, the default Automatic is used. This will set the DFT-1/2 default parameters from the table of optimized values given below.

Warning

When the DFT-1/2 method is enabled, the self-energy potential is included in the calculation of the total energy, forces and stress. However, the method is only intended as a scheme for calculating band structures, not for structural relaxations. We recommend that structural properties be calculated before applying the DFT-1/2 correction.

DFT-1/2 default parameters

The default DFT-1/2 parameters used by the Automatic keyword have been optimized against a wide range of materials, and should improve upon the standard DFT band gap in most cases. The complete list is given in the table below. As explained above, we only include elements which tend to have anionic character in compounds of interest.

Currently, only Ge (GGA) has an inner_radius different from zero, meaning that this is the only element using the shell-DFT-1/2 method as default.

Important

All elements not present in the table default to Disabled, i.e., no DFT-1/2 correction.

Important

For meta-GGA exchange-correlation functionals, all elements default to Disabled.

Table 11 Optimized DFT-1/2 parameters.

Element

LDA functionals

GGA functionals

fractional_charge

cutoff_radius | inner_radius

fractional_charge

cutoff_radius

inner_radius

C

0.3

2.5 Bohr

0.0 Bohr

0.4

2.5 Bohr

0.0 Bohr

N

0.4

3.0 Bohr

0.0 Bohr

0.4

3.0 Bohr

0.0 Bohr

O

0.5

2.5 Bohr

0.0 Bohr

0.4

3.0 Bohr

0.0 Bohr

F

0.7

3.0 Bohr

0.0 Bohr

0.6

2.5 Bohr

0.0 Bohr

Si

0.2

4.0 Bohr

0.0 Bohr

0.2

4.0 Bohr

0.0 Bohr

P

0.3

4.0 Bohr

0.0 Bohr

0.3

3.5 Bohr

0.0 Bohr

S

0.6

3.5 Bohr

0.0 Bohr

0.4

3.5 Bohr

0.0 Bohr

Cl

0.7

3.5 Bohr

0.0 Bohr

0.7

3.5 Bohr

0.0 Bohr

Zn

0.5

2.0 Bohr

0.0 Bohr

0.5

2.0 Bohr

0.0 Bohr

Ge

0.5

4.0 Bohr

0.0 Bohr

0.25

3.4 Bohr

2.35 Bohr

As

0.3

4.0 Bohr

0.0 Bohr

0.3

4.0 Bohr

0.0 Bohr

Se

0.5

3.5 Bohr

0.0 Bohr

0.5

3.5 Bohr

0.0 Bohr

Br

0.7

4.0 Bohr

0.0 Bohr

0.7

4.0 Bohr

0.0 Bohr

Sn

0.7

5.0 Bohr

0.0 Bohr

0.9

4.5 Bohr

0.0 Bohr

Sb

0.3

4.5 Bohr

0.0 Bohr

0.3

4.5 Bohr

0.0 Bohr

Te

0.4

4.0 Bohr

0.0 Bohr

0.4

4.0 Bohr

0.0 Bohr

I

0.6

4.5 Bohr

0.0 Bohr

0.6

4.0 Bohr

0.0 Bohr

Spin-orbit coupling (SOC)

Background

The spin-orbit coupling can lead to energy level splits of bandstructures and molecular energy spectra. As a rule of thumb, the magnitude of SOC scales as \(Z^{4}/n^{3}l^{2}\), where \(Z\), \(n\), and \(l\) are the atomic charge, principle, and angular quantum number of the considered energy level, respectively. The SOC is fundamentally a relativistic effect, and can be described by a second-order expansion of the relativistic Hamiltonian in terms of the fine structure constant \(\alpha\) [10].

In DFT: LCAO, the SOC is taken into account via the NormConservingPseudoPotential. For example, the SG15 pseudopotentials, shipped with QuantumATK as of the 2016 release, come in two versions for each element; a standard scalar-relativistic pseudopotential (denoted “SG15”), and one that is generated by mapping the solution to the Dirac equation, which naturally includes the SOC, to a scalar-relativistic pseudopotential (denoted “SG15-SO”):

\[V_{ps} = V_{L} + V_{NL}^{+1/2} + V_{NL}^{-1/2},\]

with a local contribution \(V_{L}\) and non-local contributions from total angular momenta \(j=l+1/2\) and \(j=l-1/2\).

Each non-local term is expanded in spin-orbit projector functions \(P_{\alpha\beta}^{l \pm 1/2, \xi}\),

\[V_{NL}^{\pm 1/2} = \sum_{l,\xi,\alpha,\beta} \nu_{l\pm 1/2,\xi} P_{\alpha\beta}^{l \pm 1/2, \xi},\]

where \(\nu_{l\pm 1/2,\xi}\) are normalization constants. The indices \(\alpha,\beta\) denote the possible spin orientations (up, down). Each nonlocal pseudopotential term requires four projector functions per expansion order \(\xi\).

Tip

Use SG15-SO or OpenMX pseudopotentials for DFT: LCAO calculations with spin-orbit.

Usage

The SOC terms are enabled by choosing a pseudopotential containing spin-orbit terms (SG15-SO or OpenMX), and appropriate settings for the ExchangeCorrelation keyword. For example, an QuantumATK Python script for silicon with SOC could contain the following specification for the calculator:

#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    BasisGGASG15SO.Silicon_Medium,
    ]
#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = SOGGA.PBE
#----------------------------------------
# Calculator
#----------------------------------------
k_point_sampling = MonkhorstPackGrid(na=9,nb=9,nc=9)
numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=k_point_sampling,
    density_mesh_cutoff=100.0*Hartree,
    )
calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )

Note the setting exchange_correlation=SOGGA.PBE, which switches on the SOC terms in the pseudopotential. Omitting this keyword or setting it to a non-SOC method, e.g. NCGGA.PBE, would still result in a calculation with the SG15-SO pseudopotential, but the SOC terms would be disregarded.

Hybrid functionals in the LCAO framework

Hybrid functionals are DFT functionals where part (or all) of the exchange in the exchange correlation functional is replaced by the exact exchange contribution. The exact exchange contribution in an LCAO basis set is defined as:

\[X_{ij} = \sum_{kl} V_{ik;jl} D_{kl}~,\]

in which \(D\) is the density matrix, and \(V\) are the 4-index coulomb integrals:

\[V_{ik;jl} = \int \int \frac{\phi_{i}^*(\mathbf{r}) \phi_{k} (\mathbf{r}) \phi_{l}^* (\mathbf{r}') \phi_{j} (\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|} d\mathbf{r} d\mathbf{r}'~.\]

In our LCAO implementation we use numerical orbitals, which makes the above integral very complicated to evaluate. To improve performance of the integral evaluations we introduce the resolution of identity (RI) approach [11] which writes a product of two orbitals as:

\[\phi_{i}^*(\mathbf{r}) \phi_{k} (\mathbf{r}) = \sum_{\mu} C_{ik}^{\mu}P_{\mu}(\mathbf{r})\]

where we have defined the auxiliary basis functions \(P_\mu\) and the auxiliary basis coefficients \(C^\mu_{ik}\). In the general RI approximation, the function \(P_\mu\) in the above expansion can be located around any of the other atoms in the system. We will limit ourselves to the case where \(\mu\) lies either on the atom \(I\) or \(K\) of the original orbitals. This approximation is called the PARI approximation [12], and is the reason for the linear computational scaling with system size.

The accuracy of the numerical evaluation of the integrals can be set by the user using the parameter class ExactExchangeParameters.

Dielectric dependent Hybrid functionals.

In HSE06 we mix 0.75 of screened PBE exchange with 0.25 of screened Fock exact exchange. The value of this optimal exchange fraction to use can be shown to be inversely proportionate to the dielectric constant, \(\epsilon\), in the material. This makes sense if we look at the two limiting cases. If we have a metal we have perfect screening and an infinite \(\epsilon\), so the exchange interaction is screened away fully and you don’t want to use Fock exchange. In the vacuum, \(\epsilon=1\), there is no screening and the correct exchange fraction is 1. The factor of 0.25 was chosen to work well for semi-conductors such as Silicon or Germanium, but is not optimal for materials with comparatively large and small band gaps.

To improve on this, various dielectric dependent hybrid functionals (DDH) have been invented. We have chosen to follow the approach suggested in [13]. In this approach they suggest to calculate an estimator for the dielectric constant:

\[\bar{g} = \frac{1}{V}\int d\mathbf{r}~\sqrt{\frac{\nabla\rho(\mathbf{r})}{\rho(\mathbf{r})}}~.\]

It is then assumed that the exchange fraction depends on this property through a fourth order polynomial:

\[\alpha = a_0 + a_1 \bar{g} + a_2 \bar{g}^2 + a_3 \bar{g}^3 + a_4 \bar{g}^4~.\]

We have fitted these parameters so they give the best possible match to the experimental band gap for a number of semi-conducting and insulating materials. You can choose this functional by setting the exchange correlation:

exchange_correlation = HybridGGA.HSE06DDH

If you want to choose a custom parameterization, fitted to a specific set you are interested in, you can do that by manually assigning the fitting parameters:

hse06_ddh_parameters = (0.16957336, 0.0, 0.0, 0.0, 0.08940146)
exchange_correlation = HybridGGA.HSE06DDH(ddh_parameters=ddh_parameters)

This functional assumes that the whole system has the same exchange fraction, so it is meaningful only when used for homogeneous bulk systems.

For interface system consisting of two or more different materials with different properties, a different approach is needed. To this end the Local DDH functional has been invented [14]. In this approach a local estimator is defined:

\[\bar{g}(\mathbf{r},\sigma) = \frac{1}{(2\pi\sigma)^{\frac{3}{2}}}\int d\mathbf{r}'~\sqrt{\frac{\nabla\rho(\mathbf{r}')}{\rho(\mathbf{r}')}}\exp\left(-\frac{|\mathbf{r}-\mathbf{r}'|^2}{2\sigma}\right)~,\]

from which an exchange fraction field is defined as:

\[a(\mathbf{r},\sigma) = a_0 + a_1 \bar{g}(\mathbf{r},\sigma) + a_2 \bar{g}(\mathbf{r},\sigma)^2 + a_3 \bar{g}(\mathbf{r},\sigma)^3 + a_4 \bar{g}(\mathbf{r},\sigma)^4.\]

This is then used to modify the 4-center coulomb integrals:

\[V_{ik;jl} = \int \int \alpha(\mathbf{r},\mathbf{r}';\sigma)\frac{\phi_{i}^*(\mathbf{r}) \phi_{k} (\mathbf{r}) \phi_{l}^* (\mathbf{r}') \phi_{j} (\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|} d\mathbf{r} d\mathbf{r}'~,\]

with

\[\alpha(\mathbf{r},\mathbf{r}';\sigma) = \sqrt{a(\mathbf{r},\sigma)a(\mathbf{r}', \sigma)}~.\]

Using this functional we can get the best results for the local density of states for interfaces, where in a \(\mathrm{Si}-\mathrm{SiO}_2\) interface a 0.25 exchange fraction will be found for the \(\mathrm{Si}\) part and a 0.35 for the \(\mathrm{SiO}_2\) part.

To use this functional you have to set the exchange correlation:

exchange_correlation = HybridGGA.HSE06LocalDDH

and once again you can set a custom parameterization by:

hse06_ddh_parameters = (0.16957336, 0.0, 0.0, 0.0, 0.08940146)
exchange_correlation = HybridGGA.HSE06LocalDDH(ddh_parameters=ddh_parameters)

A caveat for this method is that in the presence of vacuum the estimator blows up, and this has a negative influence on the convergence of the method. So we would recommend only using this functional for 3D systems.

The estimator also has difficulties with metallic systems, for which it should set the exchange fraction to zero. To improve on this we have implemented our own metallic estimator, which identifies metallic regions and sets the exchange fractions to zero in them. To use it just set the exchange correlation to:

exchange_correlation = HybridGGA.HSE06MetallicLocalDDH

When using this in a metal-semiconductor interface, it is as if you are using PBE for the metallic part and HSE06LocalDDH for the semi-conducting part.

The Auxiliary Density Matrix Method (ADMM)

Using the PARI approximation, the calculation of the exchange matrix scales linearly with system size. There is one problem, however, and that is that the scaling with basis set size is polynomial with \(\mathcal{O}(n^6)\). This means that if we need large basis sets to describe a system, the exchange calculation will become prohibitively expensive.

To overcome this issue the auxiliary density matrix method (ADMM) was developed [15], in which you run a DFT calculation with a large basis set, and use a smaller subset to compute the exchange properties.

The idea is simple, we project the density matrix of the system with the big basis set \(D\) down to a smaller basis set, and obtain the ADMM density matrix \(d\). It is then postulated that a good approximation to the exchange matrix is:

\[X_{ij}(D) \approx x_{ij}(d) + \langle \phi_i | V_x(D)| \phi_j\rangle - \langle \phi_i |V_x(d)| \phi_j \rangle~,\]

in which the latter two terms are the exchange potentials generated by the full density and the ADMM density.

Up and down projections to and from the ADMM basis are done during the diagonalization procedure. This leads to some overhead in the diagonalization procedure. In return we get a significant speedup in the exchange matrix calculation.

The accuracy of the results depends on the ADMM basis set you choose. For the PseudoDojo basis sets we have constructed custom made ADMM basis sets, which are tuned to give accurate results for band gaps and structural properties (\(\Delta\)-value) for a wide range of semi-conductors and insulators. In particular the PseudoDojo ADMM basis sets have been fitted to work best in combination with PseudoDojo Medium as main basis set, while for larger PseudoDojo basis sets we suggest to follow the hierarchy. For the FHI basis sets, you have a natural hierarchy you can choose from. For the other basis sets we don’t recommend using ADMM.

To use it, you need to set:

exact_exchange_parameters = ExactExchangeParameters(
    use_admm=True,
    basis_set_admm=BasisGGAPseudoDojo.ADMM)

numerical_accuracy_parameters = NumericalAccuracyParameters(
    exact_exchange_parameters=exact_exchange_parameters
    )

lcao_calculator = LCAOCalculator(
    exchange_correlation=HybridGGA.HSE06,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )