The ATK-DFT package


Introduction	The ATK-SE Package

Table of Contents

Introduction

Atomistix ToolKit can model the electronic properties of closed and open quantum systems based on density-functional theory (DFT) using numerical basis sets.

The key parameter in the self-consistent loop is the density matrix. For open systems, the density matrix is calculated using non-equilibrium Green's functions, while for closed or periodic systems it is calculated by diagonalizing the Kohn-Sham Hamiltonian. The Density Matrix defines the electron density, and the electron density sets up an effective potential, i.e. the Hartree and exchange-correlation potential. From the effective potential, we can obtain the Kohn-Sham Hamiltonian.

In the section called “Background information”, we describe the mathematical formalism behind the ATK-DFT Model.

	Note
	To use DFT models in ATK it is necessary to have the `ATKDFTMaster` feature in the license file (additionally, `ATKDFTSlave` features are needed to run DFT calculations in parallel). This constitutes a separate component compared to the SemiEmpirical module, and can be purchased separately.

Background information

The LCAOCalculator provides a description of the electronic structure using Density Functional Theory and norm-conserving pseudopotentials. The method is based on an expansion of the wave functions in numerical orbitals with compact support. In this section we will describe some of the the mathematical formalism behind ATK-DFT, which closely follows the implementation of Soler et al. [19].

Kohn-Sham Hamiltonian

In DFT, the electronic structure of the system is described in terms of the one-electron Kohn-Sham Hamiltonian

$\displaystyle \hat{H}_{1el} =-\frac{\hbar}{2m} \nabla^2 + V^{\rm eff}[n]({\bf r})$

In this equation the first term describes the kinetic energy of the electron. The second term describes an effective potential energy of the electron moving in the mean field from the other electrons, where the other electrons are described in terms of the total electron density $n$ .

The electron density is described in detail in the section called “Electron density” and the effective potential in the section called “Effective potential”

Solving the Kohn-Sham equations in a basis set expansion

We will now find the one-electron eigenstates of the Kohn-Sham Hamiltonian by solving the one-electron Schrödinger equation

$\displaystyle \hat{H}_{1el} \psi_{\alpha}({\bf r}) = \varepsilon_{\alpha} \psi_{\alpha}({\bf r})$

To solve this differential equation we expand the wave functions in a set of basis functions

$\displaystyle \psi_{\alpha}({\bf r})=\sum_i a_{\alpha i} \phi_i({\bf r}).$

and this modifies the differential equation into a matrix equation for determining $c_{\alpha i}$

$\displaystyle \sum_{j } H_{ij} c_{\alpha j} = \varepsilon_\alpha \sum_{j } S_{ij} c_{\alpha j}.$

The Hamiltonian matrix, $H_{ij} = \langle \phi_i | \hat{H}_{1el} | \phi_j \rangle,$ and overlap matrix $S_{ij} = \langle \phi_i | \phi_j \rangle$ are given by 3D integrals over the basis functions.

Electron density

The occupied eigenstates defines the electron density of the system

$\displaystyle n({\bf r}) = \sum_{\alpha}|\psi_\alpha({\bf r})|^2 f\left(\frac{\varepsilon_\alpha-\varepsilon_F}{k T}\right),$

where $f(x)=1/(1+e^x)$ is the Fermi function, $\varepsilon_F$ the Fermi energy, and $T$ the electron temperature.

The density is conveniently represented in terms of the density matrix

$\displaystyle n({\bf r}) = \sum_{ij} D_{ij} \phi_i({\bf r}) \phi_j({\bf r}),$

where the density matrix is given by basis set expansion coefficients

$\displaystyle D_{ij} = \sum_{\alpha}c_{\alpha i}^* c_{\alpha j} f\left(\frac{\varepsilon_\alpha-\varepsilon_F}{k T}\right)$

	Note
	For open systems (DeviceConfiguration), the density matrix is calculated using non-equilibrium Greens functions The NEGF formalism in ATK

Electron difference density

It is often convenient to compare the electron density with a superposition of atom based densities, $n^{\mathrm{atom}}({\bf r}-{\bf R}_\mu)$ , where ${\bf R}_\mu$ is the position of atom $\mu$ ,

$\displaystyle \Delta n({\bf r}) = n({\bf r}) - \sum_{\mu} n^{\mathrm{atom}}({\bf r}-{\bf R}_\mu).$

$\Delta n$ is called the electron difference density and is calculated with the ElectronDifferenceDensity keyword .

Effective potential

The Effective potential has 3 contributions:

$\displaystyle V^{\rm eff}[n] = V^{H}[n] + V^{xc}[n] + V^{ext},$

The two first terms are due to interactions with the other electrons (represented through the electron density). The first term, $V^{H}[n]$ , is the Hartree potential due to the mean-field electrostatic interaction, and the second term, $V^{xc}[n]$ , is the exchange-correlation potential caused by the quantum mechanical nature of the electrons.

$V^{ext}$ represents any other electrostatic interactions in the system. It can be separated into contributions from the ions potentials (represented by norm-conserving pseudo-potentials), and the electrostatic interaction with an external electrostatic field.

External potential

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

$\displaystyle V^{ext} = \sum_\mu V^{pseudo}_\mu + V^{gate}$

The pseudo potential has two contributions,

$\displaystyle V^{pseudo}=V^{local}+\sum_{n,n'}|\chi_n \rangle B_{n,n'}\langle \chi_{n'} |$

where the first part is a local interaction and the last part is a non-local interaction

	Note
	Due to the non-local part of the pseudo-potential, the range of the Hamiltonian is longer than twice the interaction range of the basis functions.

The $V^{gate}$ is the electrostatic potential from the external gates calculated with zero electron density. This term is calculated with the MultigridSolver and included in the ElectrostaticDifferencePotential.

Hartree potential

From the electron density, we can calculate the classical electrostatic potential, the so-called Hartree potential. The calculation of the Hartree potential is described in detail in the section called “The Hartree Potential”.

Exchange-correlation potential

In DFT, the quantum mechanical effect of the other electrons are included through the exchange-correlation energy. There exists a number of different approximations for the exchange-correlation energy and in ATK we support an extensive list of different approximations for calculating the exchange-correlation energy, see the section called “Exchange-correlation energy”

The exchange-correlation potential is the functional derivative of the exchange-correlation energy, and is a mean-field quantum mechanical interaction potential between the electrons

$\displaystyle V^{xc}[n](r) = \frac{\delta E^{XC}} {\delta n}(r)$

Total energy and forces

In DFT, the total energy is a functional of the electron density $n$ and given by

$\displaystyle E[n] = T[n] + E^{xc}[n]+E^{H}[n]+E^{ext}[n]$

where $T[n]$ is the kinetic energy of the Kohn-Sham orbitals,

$\displaystyle T[n] = \sum_{\alpha}\langle\psi_\alpha|-\frac{\hbar}{2m} \nabla^2 |\psi_\alpha\rangle f\left(\frac{\varepsilon_\alpha-\varepsilon_F}{k T}\right).$

$E^{xc}[n]$ is the exchange-correlation energy, $E^{H}[n]$ is the Hartree energy and $E^{ext}[n]$ is the interaction energy with the pseudo potential ions and other electrostatic external potentials.

The total energy is calculated with the class TotalEnergy .

The forces are calculated by differentiating the total energy with respect to the ionic coordinate of atom $i$ at position ${\bf R}_i$

$\displaystyle {\bf F}_i = -\frac{d E[n]}{d {\bf R}_i}$

Pseudopotentials

ATK uses norm-conserving pseudopotentials and is shipped with a database of pseudopotentials for the entire periodic table. This database is under review and will be updated to provide more accurate, general purpose pseudopotentials.

Through the use of the keyword NormConservingPseudoPotential it is possible to specify your own pseudo potential.

ATK uses the unified pseudo potential format (upf) defined by the PWSCF consortium. From their website, it is possible to download tools that can convert several different pseudo-potential formats into a upf file.

LCAO basis set

The electronic structure is expanded in a basis of local atomic orbitals (LCAO's)

$\displaystyle \phi_{nlm}({\bf r}) = R_{nl}(r) Y_{lm}(\hat{r}),$

where $Y_{lm}$ is a spherical harmonic and $R_{nl}$ is a radial function with compact support (i.e. exactly zero outside a confinement radius).

ATK has 3 different types of basis functions:

The basis orbitals has a number of parameters that determines 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 element. The basis sets types are listed below. The basis sets are ordered with increasing number of basis orbitals. Basis sets with more orbitals are more accurate, at the expense of increasing the computational time.

SingleZeta: One ConfinedOrbital for each occupied valence orbital in the atom.
DoubleZeta: One ConfinedOrbital and one AnalyticalSplit for each occupied valence orbital in the atom.
SingleZetaPolarized: One ConfinedOrbital for each occupied valence orbital in the atom, and one PolarizationOrbital for the first unoccupied shell in the atom.
DoubleZetaPolarized: One ConfinedOrbital and one AnalyticalSplit for each occupied valence orbital in the atom. One PolarizationOrbital for the first unoccupied shell in the atom.
DoubleZetaDoublePolarized: One ConfinedOrbital and one AnalyticalSplit for each occupied valence orbital in the atom. One PolarizationOrbital and one AnalyticalSplit for the first unoccupied shell in the atom.

The standard basis orbitals can be inspected with VNL. Open the Custom Analyzer tool . and select Analyzers → Basis Set Explorer.

DoubleZetaPolarized basis set for Silicon. The angular momentum, 14 0 bf36060e03e5988e194ab69184da668d l of the basis orbitals are evident from the asymptotic behaviour 16 0 1a68261db58957907234cb700e76e5d9 r^l at the origo. Thus the 3s basis-orbitals are red and dark-green with values 0.18 and 0.32 at origo. The 3p basis-orbitals are black and blue, and increase linearly around origo. The 4d basis-orbital is turquise and increases quadratically around origo. Note, that the basis orbitals are Gram-Schmidt orthogonalized.Gram-Schmidt orthogonalized

Figure 2: DoubleZetaPolarized basis set for Silicon. The angular momentum, $l$ of the basis orbitals are evident from the asymptotic behaviour $r^l$ at the origo. Thus the 3s basis-orbitals are red and dark-green with values 0.18 and 0.32 at origo. The 3p basis-orbitals are black and blue, and increase linearly around origo. The 4d basis-orbital is turquise and increases quadratically around origo. Note, that the basis orbitals are Gram-Schmidt orthogonalized.

Exchange-correlation energy

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

The main family of Exchange-Correlation functionals supported in ATK are Local Density Approximations(LDA) and Generalized Gradient Approximations (GGA). Each functional comes in two variants, a non-polarized and a spin-polarized version. By selecting the variant, the user will determine whatever the calculation will be spin polarized or non-polarized.

Local Density Approximation (LDA)

In LDA the exchange-correlation functional is taken as a function of the local density

$\displaystyle E^{LDA}[n] = \int n(r) \varepsilon^{LDA}(n(r)) dr,$

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

It is possible to calculate the exact exchange energy for the homogeneous electron gas, the so-called Dirac-Bloch exchange energy, which is used for all LDA functionals. The correlation energy cannot be calculated exactly, however, and a number of different approximations exists. Most of them give almost identical results; the most commonly used variant is the parametrization by Perdew and Zunger, available through the keywords LDA.PZ and LSDA.PZ.

For a list of the parametrizations for the LDA correlation available in ATK see the documentation of the ExchangeCorrelation class.

Generalized Gradient Approximation (GGA)

The Generalized Gradient Approximations (GGA) denotes local exchange-correlation functionals, where the functional depends on both the local value and the local gradient of the density

$\displaystyle E^{GGA}[n] = \int n(r) \varepsilon^{GGA}(n(r), |\Delta n(r)|) dr$

For a list of the GGA exchange-correlation functionals available in ATK see the documentation of the ExchangeCorrelation class.

Meta-GGA

Meta-GGA is a higher-level approximation to GGA which can provide a much more accurate description of the band gap of semiconductors compared to LDA and GGA. In fact, the results are often comparable to methods with a significantly higher complexity level (and associated computational cost) like GW.

A noteworthy detail about the meta-GGA exchange potential is that it can usually not be written as a derivative of an exchange functional. Therefore, meta-GGA is not recommended for geometry optimizations; for that, one should use standard GGA.

Moreover, although the method can give very accurate band gaps for semiconductors and insulators, meta-GGA typically performs rather poorly for metals. Therefore, for systems involving both metals and semiconductors/insulators, it is better to use the Hubbard U model within LDA or GGA, with the U parameters for the semiconductor elements fitted against meta-GGA results.

For a list of the meta-GGA exchange parametrizations available in ATK see the documentation of the ExchangeCorrelation class.

XC+U mean-field Hubbard term

The local approximations to the exchange-correlation energy has a number of short comings, where the most important are

Self interaction, an electron interacts with itself and this prevents electrons from localizing.
Poor approximation for excited states, band gaps are often too low.

The mean field Hubbard correction, popularly called LDA+U or XC+U, is a semi-empirical correction which tries to improve on these deficiencies of the local exchange-correlation functionals.

In the XC+U an additional energy term,

$\displaystyle E_{U} = \frac{1}{2} \sum_\mu U_\mu (n_\mu - n_\mu^2)$

is added to the Exchange-Correlation energy [26]. 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 un-occupied shell, while positive for a fractional occupied shell.

The energy is thereby lowered if states become fully occupied. This may happen if the energy levels move away from the Fermi Level, i.e. increasing the band gap, or if the broadening of the states is decreased, i.e. the electrons are localized. Thus, the Hubbard U improves on the deficiencies of the local exchange-correlation functional 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 is obtained by minimizing the total energy of the system [26].

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 implementation in ATK

The XC+U implementation in ATK comes in two variants, a standard local orbital formulation, socalled onsite representation, and the dual representation introduced by Han et al. [27]. The difference between the two implementations is in the definition of the local occupation matrix.

Onsite representation

In the onsite representation the local occupation matrix is given by

$\displaystyle n_{\mu, m m'}^\sigma = D_{\mu m, \mu m' }^\sigma$

where $D$ is the density matrix, and $\sigma$ a spin index.

Dual representation

In the dual representation the local occupation matrix is given by

$\displaystyle 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.

The onsite representation correspond to projecting into a local orbital, and since the local orbitals are non-orthogonal there are no sum rules for the occupation matrix.

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 sum up to the total number of electrons.

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

Input format for XC+U

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

The LSDAU.PZ keyword is an abbreviation of the ExchangeCorrelation keyword ExchangeCorrelation, and the specification

exchange_correlation = LSDAU.PZ

is identical to

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

Use this form with the keyword 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 a U=4.6 eV. The 3d orbitals are filled anisotropic, i.e, there are 1 electron in m=-2,-1,0,1 leaving m=2 empty. (The alternative is SphericalSymmetric which adds 4/5 electron in each orbital.)