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

$\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 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¶

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

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

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

$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 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],

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.25 3.4 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 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.

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