ATK-DFT

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).

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), while for closed or periodic systems it is calculated by diagonalization of the Kohn–Sham Hamiltonian. The electron density then sets up an effective potential, which is given by the Hartree, exchange-correlation, and external potential. Knowing the effective potential allows us to obtain the Kohn–Sham Hamiltonian.

The next section describes the mathematical formalism behind the ATK-DFT 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 ATK-DFT calculator, which closely follows the work of Soler et al. [SAG+02].

Kohn–Sham Hamiltonian

Within density funtional 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 formalism.

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 The Hartree Potential.

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

ATK uses norm-conserving pseudopotentials and is shipped with a database of pseudopotentials for the entire periodic table, see the table below. This database is reviewed and updated on a regular basis to provide increasingly accurate and general-purpose pseudopotentials.

Table 9 Pseudopotential types available for the ATK-DFT calculator.
Name Notes
FHI Trouiller–Martins type.
HGH Hartwigsen–Goedecker–Hutter type.
OMX Fully relativistic, see OpenMX database (2013).
SG15 ONCV type, scalar-relativistic (SG15) and fully relativistic (SG15-SO).

Through the use of the keyword NormConservingPseudoPotential, it is possible to specify a pseudopotential you want to use. The first three pseudopotential types are available in both LDA and GGA versions, while the SG15 pseusopotentials were generated with GGA only.

ATK 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.

Warning

Be aware that the OMX pseudopotentials, especially those including semi-core states, require a very large mesh cutoff. The required cutoff will be system-dependent, but about 400 Hartree are often appropriate.

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.

ATK comes with a number of pre-build basis sets for each chemical element and for each type of pseudopotential. The basis set types for the SG15, OMX, and FHI pseudopotentials are listed below. They are ordered with increasing number of basis orbitals; those with more orbitals are more accurate, at the expense of increasing the computational time.

SG15

Three different basis sets are available for each element when using the SG15 pseudopotentials; Medium, High, and Ultra. All three derive from the numerical atom-centered basis sets of the FHI-aims package, but have been significantly modified and optimized with respect to computational speed with the ATK-DFT calculator.

The Medium basis set is default for the SG15 pseudopotentials, and should be sufficient for most applications. However, if extreme accuracy is needed, the High and Ultra basis sets add more basis functions at the expense of increased computational load.

OMX

The OMX basis sets were developed in Refs. [Oza03] and [OK04], and are optimized towards maximum computational speed and accuracy with the OMX pseudopotentials. The ATK-DFT calculator offers three different OMX basis set sizes for most elements; Low, Medium, and High, which include increasingly more pseudo-atomic orbitals in the basis set. More details are available on the website for the OpenMX database (2013).

FHI

Four different types of basis functions can be used for the FHI pseudopotentials:

The basis sets available for FHI pseudopotentials are listed below. They are ordered with increasing number of basis orbitals – those with more orbitals are often more accurate, at the expense of increased computational load. The DZP basis set is default.

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) and the Generalized-gradient approximation (GGA). QuantumATK also supports the TB09 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 ATK-DFT 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 ATK-DFT.

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 22, Table 25, and Table 26.

Generalized-gradient approximation (GGA)

The GGA functionals is a large family of semi-local approximations for the exchange-correltion 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 22, Table 25, and Table 26.

Meta-GGA

The TB09 meta-GGA (MGGA) is a higher-level approximation that can provide a much more accurate description of semiconductor band gaps than ordinary LDA and GGA. In fact, that accuracy of semiconductor band gaps obtained with MGGA are often comparable to that of calculations with GW or hybrid functional, which have a significantly higher computational cost.

Warning

The TB09 exchange potential can not be written as a derivative of an exchange functional. First-principles forces are therefore not available, so the TB09 meta-GGA can not be used for geometry optimization; for that, one should use a standard GGA functional.

As stated above, the TB09 method can give very accurate band gaps for semiconductors and insulators, and can also be applied to systems involving both metals and semiconductors/insulators. It is even possible to tune the TB09 method such that it works optimum on both sides of an interface, see ExchangeCorrelation. Alternatively, one can use the Hubbard U model within LDA or GGA, with the U parameters for the semiconductor elements fitted against TB09 results.

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. [DBS+98] and Cococcioni et al. [CdG05], 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 [CdG05].

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. [HOY06]. 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. [FMT08], is based on the much older Slater half-occupation scheme for molecules [SJ72] (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 [Jan78]:

\[\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 Janak’s theorem 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 [FMT11].

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 [FMT08] [FMT11].

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 [FMT08] [FMT11].

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)

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.

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 10 Optimized DFT-1/2 parameters.
Element LDA functionals GGA functionals
fractional_charge cutoff_radius fractional_charge cutoff_radius
C 0.3 2.5 Bohr 0.4 2.5 Bohr
N 0.4 3.0 Bohr 0.4 3.0 Bohr
O 0.5 2.5 Bohr 0.4 3.0 Bohr
F 0.7 3.0 Bohr 0.6 2.5 Bohr
Si 0.2 4.0 Bohr 0.2 4.0 Bohr
P 0.3 4.0 Bohr 0.3 3.5 Bohr
S 0.6 3.5 Bohr 0.4 3.5 Bohr
Cl 0.7 3.5 Bohr 0.7 3.5 Bohr
Zn 0.5 2.0 Bohr 0.5 2.0 Bohr
Ge 0.5 4.0 Bohr 0.6 4.0 Bohr
As 0.3 4.0 Bohr 0.3 4.0 Bohr
Se 0.5 3.5 Bohr 0.5 3.5 Bohr
Br 0.7 4.0 Bohr 0.7 4.0 Bohr
Sn 0.7 5.0 Bohr 0.9 4.5 Bohr
Sb 0.3 4.5 Bohr 0.3 4.5 Bohr
Te 0.4 4.0 Bohr 0.4 4.0 Bohr
I 0.6 4.5 Bohr 0.6 4.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\) [FernandezSOSF06].

In ATK-DFT, 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 ATK-DFT 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.

[CdG05](1, 2) M. Cococcioni and S. de Gironcoli. Linear response approach to the calculation of the effective interaction parameters in the LDA+U method. Phys. Rev. B, 71:035105, Jan 2005. doi:10.1103/PhysRevB.71.035105.
[DBS+98]S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. Phys. Rev. B, 57:1505–1509, Jan 1998. doi:10.1103/PhysRevB.57.1505.
[FernandezSOSF06]L. Fernández-Seivane, M. A. Oliveira, S. Sanvito, and J. Ferrer. On-site approximation for spin–orbit coupling in linear combination of atomic orbitals density functional methods. J. Phys.: Condensed Matter, 18(34):7999, 2006. URL: http://stacks.iop.org/0953-8984/18/i=34/a=012.
[FMT08](1, 2, 3) Luiz G. Ferreira, Marcelo Marques, and Lara K. Teles. Approximation to density functional theory for the calculation of band gaps of semiconductors. Phys. Rev. B, 78:125116, Sep 2008. doi:10.1103/PhysRevB.78.125116.
[FMT11](1, 2, 3) Luiz G. Ferreira, Marcelo Marques, and Lara K. Teles. Slater half-occupation technique revisited: the LDA-1/2 and GGA-1/2 approaches for atomic ionization energies and band gaps in semiconductors. AIP Adv., 1(3):032119, 2011. doi:10.1063/1.3624562.
[HOY06]M. J. Han, T. Ozaki, and J. Yu. O(N) LDA+U electronic structure calculation method based on the nonorthogonal pseudoatomic orbital basis. Phys. Rev. B, 73:045110, Jan 2006. doi:10.1103/PhysRevB.73.045110.
[Jan78]J. F. Janak. Proof that $\partial E / \partial n_i = \varepsilon \,$ in density-functional theory. Phys. Rev. B, 18:7165–7168, Dec 1978. doi:10.1103/PhysRevB.18.7165.
[Oza03]T. Ozaki. Variationally optimized atomic orbitals for large-scale electronic structures. Phys. Rev. B, 67:155108, 2003. doi:10.1103/PhysRevB.67.155108.
[OK04]T. Ozaki and H. Kino. Numerical atomic basis orbitals from h to kr. Phys. Rev. B, 69:195113, 2004. doi:10.1103/PhysRevB.69.195113.
[SJ72]J. C. Slater and K. H. Johnson. Self-consistent-field $X \alpha \,$ cluster method for polyatomic molecules and solids. Phys. Rev. B, 5:844–853, Feb 1972. doi:10.1103/PhysRevB.5.844.
[SAG+02]J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal. The SIESTA method for ab initio order-N materials simulation. J. Phys.: Condensed Matter, 14(11):2745, 2002. URL: http://stacks.iop.org/0953-8984/14/i=11/a=302.