# 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.* [SAG+02].

### 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:

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,

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\):

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

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:

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:

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

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\) 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:

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:

The pseudopotential has two contributions:

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:

### Total energy and forces¶

The DFT total energy of a many-electron system is a functional of the electron density, \(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

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\):

### 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):

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,

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

#### 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,

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

#### 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

and the so-called kinetic energy density

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

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.* [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:

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

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\):

##### Dual representation¶

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

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\):

#### 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]:

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

#### Shell-DFT-1/2¶

Following the paper by Xue et al. [XYFM18] 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`

.

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\) [FernandezSOSF06].

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”):

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}\),

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:

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

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 [RRB+12] which writes a product of two orbitals as:

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 [LRW+15], 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`

.

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

[LRW+15] | Sergey V. Levchenko, Xinguo Ren, Jürgen Wieferink, Rainer Johanni, Patrick Rinke, Volker Blum, and Matthias Scheffler. Hybrid functionals for large periodic systems in an all-electron, numeric atom-centered basis framework. Computer Physics Communications, 192:60–69, 2015. URL: https://www.sciencedirect.com/science/article/pii/S001046551500079X, doi:https://doi.org/10.1016/j.cpc.2015.02.021. |

[RRB+12] | Xinguo Ren, Patrick Rinke, Volker Blum, Jürgen Wieferink, Alexandre Tkatchenko, Andrea Sanfilippo, Karsten Reuter, and Matthias Scheffler. Resolution-of-identity approach to hartree–fock, hybrid density functionals, RPA, MP2 andGWwith numeric atom-centered orbital basis functions. New Journal of Physics, 14(5):053020, may 2012. URL: https://doi.org/10.1088/1367-2630/14/5/053020, doi:10.1088/1367-2630/14/5/053020. |

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

[XYFM18] | Kan-Hao Xue, Jun-Hui Yuan, Leonardo R.C. Fonseca, and Xiang-Shui Miao. Improved lda-1/2 method for band structure calculations in covalent semiconductors. Computational Materials Science, 153:493–505, 2018. URL: https://www.sciencedirect.com/science/article/pii/S0927025618304166, doi:https://doi.org/10.1016/j.commatsci.2018.06.036. |