¹¹institutetext: Department of Physics, University of Calcutta, 92 A P C Road, Kolkata 700009, India.

Semiconductor Physics: A Density Functional Journey

Sujoy Datta Email: [email protected] Debnarayan Jana Email: [email protected]

Abstract

The journey of theoretical study on semiconductors is reviewed in a non-conventional way. We have started with the basic introduction of Hartree-Fock method and introduce the fundamentals of Density Functional Theory (DFT). From the oldest Local Density Approximations (LDA) to the most recent developments of semi-local corrections [Generalised Gradient Approximation (GGA), Meta-GGAs], hybrid functionals and orbital dependent methodologies are discussed in detail. To showcase the performance of DFT, results obtained via different approximations are compared. We indicate the success of semi-local approximations in structural properties prediction. We also show how less computationally costly but withstand architecture of some semi-local DFT methods can solve the long riddle of bandgap underestimation. In semiconductor physics, the importance of not only the band structure prediction, but also, the proper calculation of Fermi energy, and, exact finding of band alignment is argued. The comparison of Fermi energy dependent properties can channelize the theoretical studies on modern age environment-friendly researches on semiconductors, like artificial photocatalysis, energy efficient opto-electronic devices, etc. This prescription on proper choice of DFT method is potentially competent to complement the experimental findings as well as can open up a pathway of advanced semiconducting materials discoveries.
KEYWORDS: Density Functional Theory, Exchange Correlation Functional, Jacob’s Ladder, Bandgap Underestimation, Band Alignment, Lattice Constant and Bulk Modulus.

1 Introduction

The journey of the many body physics in solids has never been a smooth one. The electrons to be considered are huge in number and their mutual interactions make the situation more challenging. It is computationally impossible to solve the coupled second order differential equation exactly for any solid. The Hartree-Fock (HF) antisymmetric wavefunction, which can be exactly found for any non-interacting system or for any system under averaged interaction, can not address the correlation of electrons in a real system. The Hohenberg-Kohn (HK) theorems have provided a breathing space for theorist and the era of Density Functional Theory (DFT) has begun.

The HK theorem has provided an one-to-one correspondence between the ground state charge density and the external potential. Using this, the coupled Schrödinger equation for $N$ electrons transforms to a single particle decoupled Kohn-Sham (KS) equation. The idea of Kohn-Sham method is to find a single particle effective potential $v_{eff}({\bf r})$ such that the exact ground state density $\rho_{0}({\bf r})$ of the interacting system equals to the ground state density corresponding $v_{eff}({\bf r})$ . So, the interacting $N$ particle system is approximately transferred to a $N$ electron system $S$ under influence of averaged effective potential (S system will carry this definition throughout). Hence, the coupled equation of $N$ interacting electron system now transforms to a single particle equation.

Though this has been a game changer in material theory, however, different types of problem has appeared due to the approximations involved in forming $v_{eff}$ . Energy bandgap, which is the reason behind the vibrant properties of semiconductors and insulators is underestimated in most of the KS scheme based approximations due to inappropriate addressing of \textcolorblackexchange and correlation (xc) between electrons.

On he other hand, experimental investigations on semiconductors are prolific. It has started long ago, and even today, lots of effort is spent on producing more and more promising materials. From the introduction of DFT method, the theorists also try to work in harmony, but, in early days, the bandgap underestimation of primitive KS schemes made their research tremendously challenging. \textcolorblackAn alternative method has been proposed through the introduction of scissor operators to rigid shift the bands to match the experimental bandgaps. But this constant shift has no physical base and for those materials where no experimental bandgap data is available, this method is not useful. Furthermore, the goal of the theoretical development is not just a data matching but finding out the physical reason behind the experimental result as well as to suggest new materials with fundamental excellence and for that the solution should come from physical ground.

Perdew et al. and Sham-Schluter have provided the understanding \textcolorblackof the bandgap problem in terms of discontinuity of functional derivative of xc energy with respect to (w.r.t.) electron density, namely, the derivative discontinuity (perdew1982, ; sham1983, ). Harbola and Sahni have dug into the basic quantum mechanical theory and find out the physical interpretation of this derivative discontinuity (harbola1989, ; sahni1990, ). \textcolorblackThe next challenge has been the search for a proper approximation, which can overcome the hurdle.

A plethora of exchange-correlation functionals and/or DFT based methods have been proposed and in most of the cases, increasing accuracy demands increasing computational cost. These are arranged in a systematic way known as Jacob’s ladder (perdew2001, ) with following rungs (i) Local Density Approximation (LDA), (ii) Generalised Gradient Approximation (GGA), (iii) Meta-GGA (MGGA), (iv) Hybrid Functional, (v) Random Phase Approximation (RPA) and beyond. The first rung is local approximation, whereas, second and third rungs are semi-local ones. \textcolorblackThe higher rungs involve non-local terms.

In the set of semi-local functionals, some are semi-empirical, i.e., one or more parameter of those are found by experimental data fitting, and some are non-empirical, i.e, found by fully theoretical development. Though the structural properties of materials are almost accurately predicted by some advanced non-empirical semi-local approximations, the bandgap predicted by those are highly underestimated. \textcolorblackIt has been observed that for finite systems, the accurate treatment of xc potential, particularly in the asymptotic regions, leads to eigenvalue differences between the highest occupied (HO) and lowest unoccupied (LU) orbital energies close to true excitation energies (savin1998, ). Most of the semi-empirical corrections intended to produce exact band gap are inspired by this observation.

Among the semi-local functionals, vanLeeuwen-Baerends (vLB) (vLB, ), Becke-Johnson (BJ) (becke2006, ) and Tran-Blaha modified BJ (TBmBJ) (tran2009, ) type functionals are potential-only correction, which means that, they are not found by functional derivative of xc energy (tran2019, ).

As, HF or related schemes provide exact-exchange (EXX), so, a portion of that is mixed in various way (global, local, etc.) in hybrid functional scheme. Till date, Heyd $-$ Scuseria $-$ Ernzerhof (HSE) screened hybrid functional method is one the most successful and widely used method for bandgap calculation (HSE, ; HSE1, ). Almost exact band structure estimation using \textcolorblackthis method provides the appropriate platform for their success in calculation of optical properties of semiconductors as well.

Now, defects are ubiquitous in nature. Materials properties are highly sensitive to the formation of different types of defects. Beside the band structure, the importance of the correct determination of Fermi energy ( $E_{F}$ ) in calculation of defect formation energy is also evident (west2012, ). Exact calculation of $E_{F}$ is also a vital element in prediction of materials’ work-function. Proper prediction of work-function is the key of perfect estimation of photocatalytic activity of the semiconductors. Also, the proposal of different types of semiconductor heterojunctions and metal-semiconductor Schottky junctions (tung2014, ) rely on the exact finding of band alignment. Thus, the new era of green energy materials (clark2008, ) becomes highly dependent on the band-alignment calculation. From three and two dimensional traditional semiconductors like silicon and graphene to relatively new materials C₃N₄ , silicene, germanene, whenever the exact prediction of work-function is needed, HSE hybrid function is a natural choice for such studies and proved to match the experimental values almost exactly (west2012, ; datta2020carbon, ; singhArunima2015, ; datta2020pccp, ).

The success of different DFT methods also depends on the choice of basis set. Results differ from pseudopotential based plane-wave basis method to site-centred methods, and, sometimes the implementation of functionals is tricky. The anomaly regarding the delocalized HO and LU orbitals and very localized xc kernel in (semi-)local functionals can be solved by imposing the correct asymptotic behaviour of xc potential locally within a solid which can be implemented easily in site-centred basis-set methods. That is the situation for vLB functionals. This method which has shown success in atomic systems (banerjee1999, ) has been implemented for the bandgap calculation of solids within localized basis of atomic sphere approximation (ASA) based linearised muffin tin orbital (LMTO-ASA) package (singh2013, ; singh2017, ).

This vLB correction to LDA has found a farm footing with the introduction of self consistent full-potential N-th order muffin tin orbital (FP-NMTO) basis sets. Using this vLB-FP-NMTO method Datta et al. have showed improvement in finding the structural properties and band gaps for orthodox group IV and group III-V materials (datta2019, ). This method \textcolorblackhas further been applied on carbon nitrides (C₃N₄ ) polymorphs, and, it has helped to identify the $\gamma$ -C₃N₄ as a better candidate for photocatalysis (datta2020carbon, ).

We could start our discussion from the DFT theory, but there would be a gap in understanding the vastly used terms in DFT theory, like exchange-correlation, Fermi hole, self-interaction, and, so on. \textcolorblackDefinition of those are as much deep rooted, as the challenges of DFT approximations are. So, a discussion of Hartree-Fock theory is worthy to bridge the gap of understanding the difficulties of DFT approximate methods. This chapter starts with that, followed by a systematic discussion on methodological advancements in DFT till recent time. Then, we attempt to compare the performance of these DFT methods, talk about their success and failures, so that, finally this can provide a clear guideline on the proper choice of methods for investigation of \textcolorblackdifferent physical properties.

2 The Hartree - Fock Method

Let us start the discussion from the time independent Schrödinger equation :

\hat{H}\Phi=E\Phi~{};~{}~{}\hat{H}=\hat{T}+\hat{V}

(1)

Though the kinetic energy operator $\hat{T}$ is same for all systems, however, the potential energy operator part $\hat{V}$ of the Hamiltonian $\hat{H}$ is very much system dependent. \textcolorblackThus, the total energy $E$ becomes system dependent as well. The simplest two examples of electronic Hamiltonion are for Hydrogen atom and Helium atom expressed as:

	$\displaystyle\hat{H}_{H}$	$\displaystyle=-\frac{\hbar^{2}}{2m}\nabla^{2}-\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{r}$
	$\displaystyle\hat{H}_{He}$	$\displaystyle=-\frac{\hbar^{2}}{2m}\left(\nabla_{1}^{2}+\nabla_{2}^{2}\right)-\frac{e^{2}}{4\pi\epsilon_{0}}\left(\frac{2}{r_{1}}+\frac{2}{r_{2}}-\frac{1}{r_{12}}\right)$		(2)

Here, $m,e$ are ate electron mass and charge. \textcolorblackFor, He, there are two electrons at $r_{1}$ and $r_{2}$ w.r.t. the centre of the ion. The Hamiltonian operator for $He$ can be decoupled (in atomic unit) as:

\hat{T}=-\left(\nabla_{1}^{2}+\nabla_{2}^{2}\right)~{};~{}\hat{V}_{ie}=\hat{v}_{ie1}+\hat{v}_{ie2}=-\left(\frac{2}{r_{1}}+\frac{2}{r_{2}}\right)~{};~{}\hat{V}_{ee}=\frac{1}{r_{12}}=\frac{1}{2|{\bf r}_{1}-{\bf r}_{2}|}

(3)

The ion-electron (i-e) interaction operators ( $\hat{v}_{ie1},\hat{v}_{ie2}$ ) are operating on single electrons, and have the dimension of potential, whereas, electron-electron ( $ee$ ) interaction operator $\hat{V}_{ee}$ is a two-body energy operator, \textcolorblackwhich depends on the relative distance between the electrons. We write potential operators in small case and energy operators in capital throughout the chapter.

Due to the \textcolorblack $\hat{V}_{ee}$ term, the Schrödinger equation for two electrons can not be decoupled into two. \textcolorblackAs the atomic number increases, it becomes tougher to solve the coupled equation analytically. In addition, for molecules and condensed matter systems, there are more than one ions (atomic core) involved, so there should be ion-ion interaction as well. As the ions are heavy, they move very slowly than the electrons, so, we can begin by separating the ‘fast’ electronic from the ‘slow’ ionic degrees of freedom. This is the adiabatic or the Born-Oppenheimer approximation (born1927, ).

Situation would be simple if $\hat{V}_{ee}$ could be written in a summative form over these two electrons, so that, we could write an effective potential $\hat{v}_{s}$ replacing $\hat{v}_{ie}+\hat{v}_{ee}$ . This is an independent particle model followed by Slater-Hartree method and then $\Phi({\bf R})=\phi_{1}({\bf r}_{1})\phi_{2}({\bf r}_{2})$ , where, ${\bf R}=\{{\bf r}_{1},{\bf r}_{2}\}$ . Now, according to the Pauli exclusion principle, the wavefunction should be anti-symmetric, so, adding Fock’s contribution the Hartree-Fock wavefunction becomes: $\Phi^{HF}({\bf R})=\begin{vmatrix}\phi_{1}({\bf r}_{1})&\phi_{2}({\bf r}_{1})\\ \phi_{1}({\bf r}_{2})&\phi_{2}({\bf r}_{2})\end{vmatrix}$ . This is a $2$ electron Slater type wavefunction build up out of single electron orbitals for a system where each electron moves in an effective potential including average effect of $ee$ repulsion.

Going back to the Eq. [3], except $\hat{V}_{ee}$ , \textcolorblackall other terms can be expressed as sum over the number of electrons, i.e., they are independent-electron term. Hence, the contribution of these terms in the total energy are also additive over the electron number $N$ . The remaining e-e interaction energy term can be expressed as:

	$\displaystyle E^{HF}_{ee}=\langle\Phi^{HF}\|\hat{V}_{ee}\|\Phi^{HF}\rangle$
$\displaystyle=\iint\frac{d{\bf r}_{1}d{\bf r}_{2}}{2\|{\bf r}_{1}-{\bf r}_{2}\|}$	$\displaystyle\left[\phi^{\ast}_{1}({\bf r}_{1})\phi^{\ast}_{2}({\bf r}_{2})\{\phi_{1}({\bf r}_{1})\phi_{2}({\bf r}_{2})-\phi_{2}({\bf r}_{1})\phi_{1}({\bf r}_{2})\}\right.$
	$\displaystyle\left.-\phi^{\ast}_{2}({\bf r}_{1})\phi^{\ast}_{1}({\bf r}_{2})\{\phi_{1}({\bf r}_{1})\phi_{2}({\bf r}_{2})-\phi_{2}({\bf r}_{1})\phi_{1}({\bf r}_{2})\}\right]$
$\displaystyle=\iint\frac{d{\bf r}_{1}d{\bf r}_{2}}{2\|{\bf r}_{1}-{\bf r}_{2}\|}$	$\displaystyle\{\|\phi_{1}({\bf r}_{1})\|^{2}\|\phi_{2}({\bf r}_{2})\|^{2}+\|\phi_{2}({\bf r}_{1})\|^{2}\|\phi_{1}({\bf r}_{2})\|^{2}\}$
$\displaystyle-\iint\frac{d{\bf r}_{1}d{\bf r}_{2}}{2\|{\bf r}_{1}-{\bf r}_{2}\|}$	$\displaystyle\{\phi^{\ast}_{1}({\bf r}_{1})\phi_{1}({\bf r}_{2})\phi^{\ast}_{2}({\bf r}_{2})\phi_{2}({\bf r}_{1})-\phi^{\ast}_{1}({\bf r}_{2})\phi_{1}({\bf r}_{1})\phi^{\ast}_{2}({\bf r}_{1})\phi_{2}({\bf r}_{2})\}$	(4)

See, $\phi_{1}$ at ${\bf r}_{1}$ and ${\bf r}_{2}$ , i.e., $\phi_{1}({\bf r}_{1})$ and $\phi_{1}({\bf r}_{2})$ are same electronic state at two \textcolorblackpositions, so, $\phi_{1}({\bf r}_{1})\phi_{1}({\bf r}_{2})$ type term represents $e-e$ self-interaction, which is unphysical. In the terms of Eq. [2], self-interaction term, i.e., $|\phi_{1}({\bf r}_{1})|^{2}|\phi_{1}({\bf r}_{2})|^{2}+|\phi_{2}({\bf r}_{1})|^{2}|\phi_{2}({\bf r}_{2})|^{2}$ is absent, so, HF theory is a theory of exact-exchange and is self-interaction corrected (SIC) from the very beginning.

2.1 Hartree - Fock Method for $N$ electron system

Let us now generalize the idea to electronic many-body systems beyond He atom. \textcolorblackThe general time-independent Schrödinger equation involving many-electron Hamiltonian is:

	$\displaystyle\hat{H}\Psi({\bf X})$	$\displaystyle=E\Psi({\bf X})\text{ ; }{\bf X}=\{{\bf r}_{1}\bm{\sigma}_{1},{\bf r}_{2}\bm{\sigma}_{2},..,{\bf r}_{N}\bm{\sigma}_{N}\}~{}~{}\&~{}~{}~{}\hat{H}=\hat{T}+\hat{V}_{ext}+\hat{V}_{ee}$		(5)
	where,	$\displaystyle\hat{T}=\sum\limits_{i=1}^{N}\left(-\frac{1}{2}\nabla_{i}^{2}\right)\text{ ; }\hat{V}_{ext}=\sum\limits_{i=1}^{N}{v_{ext}({\bf r}_{i})}\text{ ; }\hat{V}_{ee}=\frac{1}{2}\sum\limits_{\begin{subarray}{c}i,j=1\\ i\neq j\end{subarray}}^{N}\frac{1}{r_{ij}}$		(6)

This $v_{ext}({\bf r}_{i})$ is the external potential which includes the ion-electron ( $ie$ ) interaction, \textcolorblackand, local in nature. ¹¹1An operator $\hat{A}$ is local if it follows $\hat{A}({\bf r},{\bf r}^{\prime})=\hat{A}({\bf r}^{\prime})\delta({\bf r}-{\bf r}^{\prime})$ . Simply, it means that, the operator acting at any point ${\bf r}$ does not depend on anything at ${\bf r}^{\prime}$ . Potential part of one-body Hamiltonian is often local. The many-electron wave-function is a function of combination of coordinates ( ${\bf R}=\{{\bf r}_{i}\}=\{{\bf r}_{1},{\bf r}_{2},...,{\bf r}_{N}\}$ ) and spins ( $\sigma$ $=\{\sigma_{i}\}=\{\sigma_{1},\sigma_{2},...,\sigma_{N}\}$ ) of individual electrons. The normalized spin-orbital wavefunction of HF system $\Psi^{HF}({\bf X})$ is of single Slater determinant ( $\Psi^{S}$ ) form:

\displaystyle\Psi^{HF}({\bf X})\equiv\Psi^{S}({\bf X})=\frac{1}{\sqrt{N!}}\begin{vmatrix}\psi_{1}({\bf x}_{1})&\psi_{2}({\bf x}_{1})&\cdots&\psi_{N}({\bf x}_{1})\\ \psi_{1}({\bf x}_{2})&\psi_{2}({\bf x}_{2})&\cdots&\psi_{N}({\bf x}_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \psi_{1}({\bf x}_{N})&\psi_{2}({\bf x}_{N})&\cdots&\psi_{N}({\bf x}_{N})\end{vmatrix};~{}{\bf x_{n}}={{\bf r}_{n}\sigma_{n}}

(7)

Here, $\frac{1}{\sqrt{N!}}$ is the normalisation constant and the spin-orbitals follow orthonormality condition: $\int\psi^{\ast}_{i}({\bf x}_{n})\psi_{j}({\bf x}_{n})d{\bf x}_{n}=\delta_{ij}$ . Since spin-orbitals corresponding spin up states are automatically orthonormal to those of spin down states, so, this orthonormality condition reduces to the orthonormality of space-only orbitals ( $\{\phi_{i}({\bf r}_{n})\}$ ) having parallel spin. Writing the spin-orbitals $\psi_{i}({\bf x}_{n})=\phi_{i}({\bf r}_{n})\chi_{i}(\sigma_{n})$ and as the spin functions are orthonormal, we can build the spatial only HF wavefunction $\Phi^{HF}({\bf R})$ through replacing $\{\psi_{j}({\bf x}_{n})\}$ by $\{\phi_{j}({\bf r}_{n})\}$ .

Note that, the Slater type wavefunction requires separation of variable for each spin-orbitals which is only possible if there is no two-body term in the Hamiltonian. Two options are there, either $\hat{V}_{ee}=0$ , or, $\hat{V}_{ee}$ is approximated as an averaged potential. In HF scheme, we try to build a Slater type wavefunction for system having $\hat{V}_{ee}\neq 0$ . Using this trial wavefunction $\Psi^{HF}({\bf X})\equiv\Psi^{S}({\bf X})$ , let us express the e-e interaction energy $E_{ee}$ (hartree1928, ; parr1980, ).

	$\displaystyle E_{ee}^{HF}=\langle\Psi^{HF}\|\hat{V}_{ee}\|\Psi^{HF}\rangle=\frac{1}{2}$	$\displaystyle\sum_{i,j=1}^{N}\left(J_{ij}-K_{ij}\right)=E_{H}+E_{x}$		(8)
	$\displaystyle J_{ij}=\iint\frac{\psi_{i}^{\ast}({\bf x})\psi_{i}({\bf x})\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{j}({\bf x}^{\prime})}{\|{\bf r}-{\bf r}^{\prime}\|}d{\bf x}d{\bf x}^{\prime}$	$\displaystyle;K_{ij}=\iint\frac{\psi_{i}^{\ast}({\bf x})\psi_{i}({\bf x}^{\prime})\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{j}({\bf x})}{\|{\bf r}-{\bf r}^{\prime}\|}d{\bf x}d{\bf x}^{\prime}$

Here, $J_{ij}$ and $K_{ij}$ are Coulomb interaction integral and exchange energy integral and the sums are over all $\{i,j\}$ , as the self terms: $J_{ii}=K_{ii}$ cancel each other.

The Hartree energy $(E_{H})$ and corresponding Hartree potential $v_{H}({\bf r})$ are expressed by summing over $J_{ij}$ through the definition of density $\rho({\bf r})=\sum\limits_{i=1}^{N}\psi^{\ast}_{i}({\bf x})\psi_{i}({\bf x})$ as:

\displaystyle E_{H}=\frac{1}{2}\sum_{i,j=1}^{N}J_{ij}=\frac{1}{2}\iint\frac{\rho({\bf r})\rho({\bf r}^{\prime})}{|{\bf r}-{\bf r}^{\prime}|}d{\bf r}d{\bf r}^{\prime}~{};~{}~{}v_{H}({\bf r})=\int\frac{\rho({\bf r}^{\prime})}{|{\bf r}-{\bf r}^{\prime}|}d{\bf r}^{\prime}

(9)

{svgraybox}

B1: Electron Density and Density Matrix

The number of electrons per unit volume at a position ${\bf r}_{1}$ in a given state is defined as the electron density for that state at ${\bf r}_{1}$ :

\rho({\bf r}_{1})=N\idotsint\Psi^{\ast}({\bf X})\Psi({\bf X})d\sigma_{1}d{\bf x}_{2}\cdots d{\bf x}_{N}~{};~{}\int\rho({\bf r}_{1})d{\bf r}_{1}=N

(10)

For Slater type wavefunction $\Psi^{S}({\bf X})$ (as in HF scheme)

\displaystyle\int[\cdots]d{\bf x}\equiv\sum_{\sigma}\int[\cdots]d{\bf r}\Rightarrow\rho({\bf r})=\sum_{\sigma}\sum_{i}^{occ}|\psi_{i}({\bf r}_{i}\sigma)|^{2}

(11)

The reduced density matrices of first and second order are:

		$\displaystyle\gamma_{1}({\bf x}_{1},{\bf x}_{1}^{\prime})=N\int\Psi^{\ast}({\bf x}_{1},{\bf x}_{2},{\bf x}_{3},\cdots,{\bf x}_{N})\Psi({\bf x}_{1}^{\prime},{\bf x}_{2},{\bf x}_{3},\cdots,{\bf x}_{N})d{\bf x}_{2}\cdots d{\bf x}_{N}$		(12)
		$\displaystyle\gamma_{2}({\bf x}_{1}{\bf x}_{2},{\bf x}_{1}^{\prime}{\bf x}_{2}^{\prime})={N\choose 2}\int\Psi^{\ast}({\bf x}_{1},{\bf x}_{2},{\bf x}_{3},\cdots,{\bf x}_{N})\Psi({\bf x}_{1}^{\prime},{\bf x}_{2}^{\prime},{\bf x}_{3},\cdots,{\bf x}_{N})d{\bf x}_{3}\cdots d{\bf x}_{N}$		(12)

The first-order and second-order spinless reduced density matrices are:

	$\displaystyle\rho_{1}({\bf r}_{1},{\bf r}_{1}^{\prime})$	$\displaystyle=\int\gamma_{1}({\bf r}_{1}\sigma_{1},{\bf r}_{1}^{\prime}\sigma_{1})d\sigma_{1}\Rightarrow\rho({\bf r}_{1})=\rho_{1}({\bf r}_{1},{\bf r}_{1})$		(13)
	$\displaystyle\rho_{2}({\bf r}_{1}{\bf r}_{2},{\bf r}_{1}^{\prime}{\bf r}_{2}^{\prime})$	$\displaystyle=\iint\gamma_{2}({\bf r}_{1}\sigma_{1}{\bf r}_{2}\sigma_{2},{\bf r}_{1}^{\prime}\sigma_{1}{\bf r}_{2}^{\prime}\sigma_{2})d\sigma_{1}d\sigma_{2}$		(13)

\text{For Slater wavefn.}~{}\rho_{1}({\bf r}_{1},{\bf r}_{1}^{\prime})=\sum_{\sigma_{1}}\sum_{i}^{occ}\psi_{i}^{\ast}({\bf r}_{1}\sigma_{1})\psi_{i}({\bf r}_{1}^{\prime}\sigma_{1})

(14)

Here, ${\bf r}_{1},{\bf r}_{1}^{\prime},\sigma_{1}$ are dummy index, so, can simply be replaced by ${\bf r},{\bf r}^{\prime},\sigma$ . $\sum\limits_{\sigma_{1}}$ is for up and down spins, $\sum\limits_{i}^{occ}$ is for all occupied orbitals and both sum yields total $N$ number of sum over $N$ spin-orbitals. For closed shell, every orbital is doubly spin degenerate. As a result, the summation over $i$ runs from $1$ to $N/2$ occupies orbitals. Another thing to notice is that, in spinless reduced density matrices, integrations are done over all spins, $\sigma_{2}\cdots\sigma_{N}$ in Eq. [12] and on $\sigma_{1}$ in Eq. [13], so, we write in term of only space orbitals as:

\rho({\bf r})=2\sum_{i}^{N/2}|\phi_{i}({\bf r})|^{2}=\rho_{\uparrow}({\bf r})+\rho_{\downarrow}({\bf r})~{}~{}\&~{}~{}\rho_{1}({\bf r},{\bf r}^{\prime})=2\sum_{i}^{N/2}\phi_{i}^{\ast}({\bf r})\phi_{i}({\bf r}^{\prime})

(15)

2.2 Hartree - Fock Equation

If we want to use the HF scheme for the interacting electronic system, then, we have to start with taking a trial wavefunction $\Psi({\bf X})$ of Slater determinant type (Eq. [7]). Using the ortho-normalization condition at any position ${\bf r}$ , $\int\psi^{\ast}_{i}({\bf x})\psi_{j}({\bf x})d{\bf x}=\delta_{ij}$ and $\frac{\partial\langle\hat{H}\rangle}{\partial\psi^{\ast}_{i}({\bf x})}=0$ , we reach to HF equation applying Lagrange’s multiplier method:

	$\displaystyle\left(-\frac{1}{2}\nabla_{i}^{2}\right)\psi_{i}({\bf x})+v_{ext}({\bf r})\psi_{i}({\bf x})$	$\displaystyle+\left[\sum\limits_{j=1}^{N}\int\frac{\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{j}({\bf x}^{\prime})}{\|{\bf r}-{{\bf r}^{\prime}}\|}d{\bf x}^{\prime}\right]\psi_{i}({\bf x})$
		$\displaystyle-\left[\sum\limits_{j=1}^{N}\int\frac{\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{i}({\bf x}^{\prime})}{\|{\bf r}-{{\bf r}^{\prime}}\|}d{\bf x}^{\prime}\right]\psi_{j}({\bf x})=\epsilon_{i}\psi_{i}({\bf x})$		(16)

For a closed shell system having even number of electrons the $N$ spin-orbitals can be divided into $N/2$ space-only orbitals, and for each of these, two spin (up and down) states.

\textcolor

blackThe Kinetic energy contribution to the total energy can be calculated exactly as:

\displaystyle T=\langle\Psi^{HF}|\hat{T}|\Psi^{HF}\rangle=-\frac{1}{2}\sum_{\sigma}\sum_{i}\int\psi^{\ast}_{i}({\bf r},\sigma)\nabla_{i}^{2}\psi_{i}({\bf r},\sigma)d{\bf r}

(17)

If we take inner product of HF Eq. [2.2] with $\psi_{i}({\bf x})$ and sum over $\{i\}$ then the HF total energy $E[\Psi^{HF}]$ in terms of orbital energy $\epsilon_{i}$ is found using Eq. [8] as:

	$\displaystyle\langle\Psi^{HF}\|\hat{T}+\hat{V}_{ext}\|\Psi^{HF}\rangle$	$\displaystyle+2E_{ee}^{HF}=\sum_{i}^{N}\epsilon_{i}$
	$\displaystyle\Rightarrow E[\Psi^{HF}]=\langle\Psi^{HF}\|\hat{T}$	$\displaystyle+\hat{V}_{ext}\|\Psi^{HF}\rangle+E_{ee}^{HF}=\sum_{i}^{N}\epsilon_{i}-E_{ee}^{HF}$		(18)

From the last term of Eq. [2.2] the orbital dependent exchange operator $\hat{v}_{x,i}({\bf r})$ and orbital dependent Fermi hole $\rho_{x,i}({\bf r},{{\bf r}^{\prime}})$ are defined as(harbola1989, ; sahni1990, ):

	$\displaystyle\hat{v}_{x,i}({\bf x})\psi_{i}({\bf x})=-\left[\sum\limits_{j=1}^{N}\int\frac{\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{i}({\bf x}^{\prime})\psi_{j}({\bf x})}{\psi_{i}({\bf x})\|{\bf r}-{{\bf r}^{\prime}}\|}d{\bf x}^{\prime}\right]\psi_{i}({\bf x})=\left[\int\frac{\rho_{x,i}({\bf r},{{\bf r}^{\prime}})}{\|{\bf r}-{{\bf r}^{\prime}}\|}d{\bf r}^{\prime}\right]\psi_{i}({\bf x})$		(19)
	$\displaystyle\text{Using}~{}~{}\int d{\bf x}\equiv\sum_{\sigma}\int d{\bf r}~{};~{}~{}{\bf x}={\bf r}\sigma,{\bf x}^{\prime}={\bf r}^{\prime}\sigma~{}~{}\text{we find:}$
	$\displaystyle\rho_{x,i}({\bf r},{{\bf r}^{\prime}})=-\frac{\sum\limits_{\sigma}\sum\limits_{j}^{occ}\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{i}({\bf x}^{\prime})\psi_{j}({\bf x})}{\psi_{i}({\bf x})}=-\sum_{\sigma}\sum_{j}\frac{\psi_{i}^{\ast}({\bf x})\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{i}({\bf x}^{\prime})\psi_{j}({\bf x})}{\psi_{i}^{\ast}({\bf x})\psi_{i}({\bf x})}$		(20)
	$\displaystyle~{}\text{ and, }~{}\int\rho_{x,i}({\bf r},{{\bf r}^{\prime}})d{{\bf r}^{\prime}}=-1$		(21)

The double sum yields $N$ electronic spin-orbitals with sum over $\{\sigma\}$ running over up $+$ down spins and sum over $\{i\}$ over all occupied orbitals corresponding each type of spin. \textcolorblackEq. [19] shows that the exchange potential operator $\hat{v}_{x,i}$ is not multiplicative, i.e., non-local operator. Furthermore, due to the dynamical, i.e., orbital dependent nature of hole, charge distribution $\rho_{x,i}$ , and, the corresponding potentials are orbital dependent. Therefore, HF equation and theory is orbital dependent, and, the calculation becomes rigorous.

2.3 Hartree - Fock - Slater Equation

Till now, \textcolorblackwe have not talked about any approximation within the HF theory. Let us note the condition, $\int\rho_{x,i}({\bf r},{{\bf r}^{\prime}})d{{\bf r}^{\prime}}=-1$ is true for each electron (occupying an spin-orbital $\psi_{i}$ ) at position ${\bf r}$ . So, an weighted averaging of $\rho_{x,i}({\bf r},{{\bf r}^{\prime}})$ over all orbitals can provide the expression of the orbital-independent, or, simple Fermi hole $\rho_{x}$ .

	$\displaystyle\rho_{x}({\bf r},{\bf r}^{\prime})=\sum_{\sigma}\sum_{i}^{occ}\rho_{x,i}({\bf r},{{\bf r}^{\prime}})p_{i}({\bf r})=-\left[\frac{\sum_{\sigma}\sum_{i,j}^{occ}\psi_{i}^{\ast}({\bf x})\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{i}({\bf x}^{\prime})\psi_{j}({\bf x})}{\sum_{\sigma}\sum_{k}\psi^{\ast}_{k}({\bf x})\psi_{k}({\bf x})}\right]$		(22)
	$\displaystyle\text{where,}~{}p_{i}({\bf r})=\frac{\psi^{\ast}_{i}({\bf x})\psi_{i}({\bf x})}{\sum\limits_{\sigma}\sum\limits_{k}^{occ}\psi^{\ast}_{k}({\bf x})\psi_{k}({\bf x})}~{}\text{is the probability}~{}\&~{}\int\rho_{x}({\bf r},{{\bf r}^{\prime}})d{{\bf r}^{\prime}}=-1$		(23)

Now, for even $N$ numbered electronic system, forming closed shells, we can define the Fermi hole using the spinless reduced density matrix using Eq. [14] and [15] as:

\displaystyle\rho_{x}({\bf r},{\bf r}^{\prime})=-\frac{|\rho_{1}({\bf r},{\bf r}^{\prime})|^{2}/4}{\rho({\bf r})/2}=-\frac{|\rho_{1}({\bf r},{\bf r}^{\prime})|^{2}}{2\rho({\bf r})}

[closed shell]

(24)

If we approximate $\rho_{x,i}({\bf r},{{\bf r}^{\prime}})$ of Eq. [19] by $\rho_{x}({\bf r},{{\bf r}^{\prime}})$ then, the exchange potential operator $\hat{v}_{x,i}({\bf r})$ becomes multiplicative or local operator $v_{x}^{S}({\bf r})$ and the Hartree-Fock-Slater equation from Eq. [2.2] using Eq. [9] becomes:

\left[-\frac{1}{2}\nabla_{i}^{2}+v_{ext}({\bf r})+v_{H}({\bf r})+v_{x}^{S}({\bf r})\right]\psi_{i}({\bf x})=\epsilon_{i}\psi_{i}({\bf x})

(25)

Where, the Slater potential $v^{S}_{x}({\bf r})$ is defined in term of Fermi hole and exchange only pair-correlation function $g_{x}({\bf r},{\bf r}^{\prime})$ as:

v_{x}^{S}=\int\frac{\rho_{x}({\bf r},{{\bf r}^{\prime}})}{|{\bf r}-{{\bf r}^{\prime}}|}d{\bf r}^{\prime}=\int\frac{\rho({\bf r}^{\prime})\{g_{x}({\bf r},{{\bf r}^{\prime}})-1\}}{|{\bf r}-{{\bf r}^{\prime}}|}d{\bf r}^{\prime}~{};~{}g_{x}({\bf r},{\bf r}^{\prime})=1+\frac{\rho_{x}({\bf r},{\bf r}^{\prime})}{\rho({\bf r}^{\prime})}

(26)

Using $\rho_{x}({\bf r},{{\bf r}^{\prime}})$ , we can also define the Pauli exchange energy $E_{x}^{HF}$ using Eq. [8]:

E_{x}^{HF}=-\frac{1}{2}\sum_{i,j=1}^{N}K_{ij}=\frac{1}{2}\iint\frac{\rho({\bf r})\rho_{x}({\bf r},{{\bf r}^{\prime}})}{|{\bf r}-{\bf r}^{\prime}|}d{\bf r}d{\bf r}^{\prime}

(27)

This is how a many electron Schrödinger equation is developed in Slater-Hartree-Fock theory without any two-body term in Hamiltonian. As, electronic many body wavefunction acting under averaged effective potential can be written as Slater determinant, starting the calculation with $\Psi^{S}$ is justified. As for the system, the Fermi hole is developed due to the self interaction and Pauli exchange effect, the exchange-only interpretation of Fermi hole $\rho_{x}$ is justified as well.

An interesting system is a system of (spin polarised) single electron. For that no sum is needed in Eq. [22] and $i,j,k,\sigma=1$ in the definition of Fermi hole and density. So, for this special case $\left.\rho_{x}({\bf r},{{\bf r}^{\prime}})\right|_{1}=-\left.\rho({\bf r}^{\prime})\right|_{1}$ and the exchange energy $\left.E_{x}^{HF}\right|_{1}=-\left.E_{H}\right|_{1}\Rightarrow\left.E_{ee}^{HF}\right|_{1}=0$ , which implies:

Hartree-Fock theory is single electron self-interaction free.

NOTE: If the wavefunction $\Psi({\bf X})$ is found for interacting system, i.e, as a solution of $\hat{H}$ including $\hat{V}_{ee}$ , then the coulomb correlation effect also counts in. For such interacting electronic system, pair-correlation function $g_{xc}({\bf r},{\bf r}^{\prime})$ and Fermi-Coulomb hole $\rho_{xc}({\bf r},{\bf r}^{\prime})$ includes Coulomb interaction between electrons.

$\displaystyle\rho({\bf r}^{\prime})~{}g_{xc}({\bf r},{\bf r}^{\prime})$	$\displaystyle=\rho({\bf r}^{\prime})+\rho_{xc}({\bf r},{\bf r}^{\prime})~{};~{}~{}\int\rho_{xc}({\bf r},{\bf r}^{\prime})d{\bf r}^{\prime}=-1$	(28)
$\displaystyle\therefore E_{ee}$	$\displaystyle=\frac{1}{2}\iint\frac{\rho({\bf r})\rho({\bf r}^{\prime})g_{xc}({\bf r},{\bf r}^{\prime})}{\|{\bf r}-{\bf r}^{\prime}\|}d{\bf r}d{\bf r}^{\prime}$
	$\displaystyle=\underbrace{\frac{1}{2}\iint\frac{\rho({\bf r})\rho({\bf r}^{\prime})}{\|{\bf r}-{\bf r}^{\prime}\|}d{\bf r}d{\bf r}^{\prime}}_{\begin{subarray}{c}\textbf{Hartree energy ($E_{H}$)}\end{subarray}}+\underbrace{\frac{1}{2}\iint\frac{\rho({\bf r})\rho_{xc}({\bf r},{\bf r}^{\prime})}{\|{\bf r}-{\bf r}^{\prime}\|}d{\bf r}d{\bf r}^{\prime}}_{\begin{subarray}{c}\textbf{Pauli-Coulomb energy ($E_{xc}$)}\end{subarray}}$	(29)

3 Hohenberg-Kohn (HK) Theorem

•

First theorem: The non-degenerate ground state density $\rho_{0}({\bf r})$ determines the external potential $v_{ext}({\bf r})$ to within a trivial additive constant.
•

Second Theorem: The non-degenerate ground state density $\rho_{0}({\bf r})$ can be determined from the ground state energy functional $E[\rho_{0}]$ via the variational principle by varying only the density.

The ground state non-degenerate density $\rho_{0}({\bf r})$ exactly determines the electron number $N$ , the ground state wave-function $\Psi[\rho_{0}({\bf r})]$ and all other electronic properties. \textcolorblackNot necessarily, $\Psi$ is a single Slater type wavefunction $\Psi^{S}$ of Eq. [7], but, can be any general $N$ electron wavefunction. In addition according to the first HK theorem, $\rho_{0}({\bf r})$ determines $v_{ext}({\bf r})$ as well (hohenberg1964, ; kohn1965, ). So, what does it mean? We can solve $N$ -electron Schrödinger equation having $\hat{V}_{ext}({\bf r})$ in Hamiltonian to find the ground state wavefunction $\Psi$ (map $A$ : $\hat{V}_{ext}({\bf r})\rightarrow\Psi$ ), and then, use $\Psi$ to find the ground state density $\rho_{0}({\bf r})$ (map $B$ : $\Psi\rightarrow\rho_{0}({\bf r})$ ). This $AB$ mapping is what the first theorem provides. It also provides the definition of $v$ -representability. ²²2If we solve Schrödinger equation with any $v_{ext}({\bf r})$ to find the antisymmetric ground state wave function $\Psi$ and using this we find $\rho_{0}({\bf r})$ , then $\rho_{0}({\bf r})$ is v-representable.

Now the first question is the inverse mapping $(AB)^{-1}$ exists? The answer is yes, and, the existence can be proved easily (parr1980, ; sahni2016, ). Second question is how this can be found? The $A^{-1}$ can be found using the inverse of the Schrödinger equation to within a constant. Another way is to use quantal Newtonian first law using force filed. This is where the Quantal DFT (Q-DFT) starts to diverge, and detail study can be found in (sahni2016, ). As, there can be infinite number of antisymmetric wavefunctions that integrate to form the ground state density $\rho_{0}({\bf r})$ , the inverse mapping $B^{-1}$ is not that straightforward as $A^{-1}$ , and can be found using constrained search approach (levy1985, ).

Let, $\rho^{t}_{0}({\bf r})$ is the trial ground state density. According to first theorem, this $\rho^{t}_{0}({\bf r})$ exactly determines the trial wave-function $\Psi^{t}[\rho^{t}_{0}]$ and external potential $v^{t}_{ext}({\bf r})$ .

\text{Now,}~{}~{}E^{t}=E[\rho^{t}_{0}]=\langle\Psi^{t}|\hat{H}|\Psi^{t}\rangle\geqslant E\text{ (ground state energy)}

(30)

As, no other variable is involved in the process of this equation except the density, so, variational principle can be applied using $\int\rho_{0}({\bf r})d{\bf r}=N$ as:

\delta\{E[\rho_{0}]-\mu(\int\rho_{0}({\bf r})d{\bf r}-N)\}=0~{}~{}~{}~{}~{}~{}~{}\Rightarrow\frac{\delta E[\rho_{0}]}{\delta\rho_{0}({\bf r})}=\mu

(31)

This is the Euler-Lagrange equation. Let, $\rho^{(N)}_{0}({\bf r})$ is the density (solution of Eq. [31]) for $N$ electron system with ground state energy $E[\rho^{(N)}_{0}]$ . If the charge is changed as $N\rightarrow N+f$ , then the energy difference is:

$\displaystyle E^{(N+f)}$	$\displaystyle-E^{(N)}=\int\left.\frac{\delta E[\rho]}{\delta\rho({\bf r})}\right\|_{\rho^{(N)}_{0}}\{\rho^{(N+f)}_{0}({\bf r})-\rho^{(N)}_{0}({\bf r})\}d{\bf r}$
	$\displaystyle=\mu(N)\int\{\rho^{(N+f)}_{0}({\bf r})-\rho^{(N)}_{0}({\bf r})\}d{\bf r}=\mu(N)\{(N+f)-N\}=f\phantom{5}\mu(N)$	(32)
	$\displaystyle\text{For }f\to 0:\lim_{f\to 0}\frac{E^{(N+f)}-E^{(N)}}{(N+f)-N}=\frac{\partial E[N]}{\partial N}=\mu^{(N)}$	(33)

So, $\mu^{(N)}$ is the energy to change the electron number, so, $\mu^{(N)}$ represents the chemical potential for $N$ electron system.

The ionization energy ( $I$ ) for removing one electron from $N$ electron system and electron affinity ( $A$ ) for adding one electron in that are, respectively:

\mu^{(N)}_{-1}=\frac{E^{(N-1)}-E^{(N)}}{(N-1)-N}=-I^{(N)}~{}~{}~{}\&~{}~{}~{}\mu^{(N)}_{+1}=\frac{E^{(N+1)}-E^{(N)}}{(N+1)-N}=-A^{(N)}

(34)

\textcolor

blackConsequently, the fundamental (Optical) band gap $E_{g}$ for isolated system of $N$ electron is:

E_{g}=I^{(N)}-A^{(N)}=\mu^{(N)}_{+1}-\mu^{(N)}_{-1}

(35)

In the Hamiltonian operator, $\hat{T}$ and $\hat{V}_{ee}$ have same form for every many electronic system. The external potential energy operator $\hat{V}_{ext}=\hat{V}_{ext}[\rho]$ is system dependent, as it depends on the external potential, i.e., knowledge of nuclei or other sources have to be known. So, we can decompose the total energy as:

	$\displaystyle E[\rho_{0}]=\langle\Psi[\rho_{0}]\|\hat{H}\|\Psi[\rho_{0}]\rangle=$	$\displaystyle F^{HK}[\rho_{0}]+\int\rho_{0}({\bf r})v_{ext}({\bf r})d{\bf r}$		(36)
	$\displaystyle\text{where,}~{}~{}F^{HK}[\rho_{0}]=$	$\displaystyle\langle\Psi[\rho_{0}]\|(\hat{T}+\hat{V}_{ee})\|\Psi[\rho_{0}]\rangle$		(37)

This $F^{HK}[\rho_{0}]$ is the universal functional of density, which is same for all electronic system. Using Eq. [36], the Euler-Lagrange Equation [31] becomes :

\frac{\delta F^{HK}[\rho_{0}]}{\delta\rho_{0}({\bf r})}+v_{ext}({\bf r})=\mu

(38)

{svgraybox}

B2: Fractional Electron Number and Derivative Discontinuity

The derivation of Eq. [33] mathematically uses a condition of fractional electron number $N$ . In an open system where electron exchange is possible between atoms, the time average of charge can be fractional. This has been shown by Perdew, Parr, Levy and Balduz (PPLB) (perdew1982, ). The extension to fractional spins and the further combination of fractional charges and spins have been done later (cohen2008, ). As quantum mechanical systems with fluctuating number of particles are described by ensemble mixture, using Levy’s constrained search formalism (levy1985, ) the charge density, energy and chemical potentials are found as:

\rho_{0}^{(N+f)}=(1-f)\rho_{0}^{(N)}+f\rho_{0}^{(N+1)}~{}~{}\&~{}~{}~{}E^{(N+f)}=(1-f)E^{(N)}+fE^{(N+1)}

(39)

	$\displaystyle\therefore\mu^{(N)}_{+}$	$\displaystyle=\left.\frac{\partial E[N]}{\partial N}\right\|_{N+}=E^{(N+1)}-E^{(N)}=-A^{(N)}\text{ for }0<f\leq+1$		(40)
	$\displaystyle\mu^{(N)}_{-}$	$\displaystyle=\left.\frac{\partial E[N]}{\partial N}\right\|_{N-}=E^{(N)}-E^{(N-1)}=-I^{(N)}\text{ for }-1\leq f<0$		(40)

\Delta\mu^{(N)}=I^{(N)}-A^{(N)}=E^{(N+1)}+E^{(N-1)}-2E^{(N)}

(41)

Hence, $\mu^{(N)}$ jumps discontinuously by an amount $\Delta\mu^{(N)}$ whenever crossing the integer electron number $N$ . The functional derivative of total energy $\left(\frac{\partial E[N]}{\partial N}\right)$ \textcolorblackshows similar nature. Hence, the derivative discontinuity occurs at each integer electron number. So, we see that the famous derivative discontinuity is a property of all methods following HK theorem, hence, approximations in Kohn-Sham DFT should satisfy that property. From Eq. [39], the piecewise linearity condition is found as $E^{(N+f)}=E^{(N)}-f\{E^{(N)}-E^{(N+1)}\}$ which further implies that $E^{(N+f)}+E^{(N-f)}>2E^{(N)}\Rightarrow I^{(N+1)}<I^{(N)}$ . Fig.1 may provide more insight.

Refer to caption — Figure 1: Variation of energy $E^{(N+f)}$ with fractional electron number $N$ showing the piecewise linearity. $I^{(N-1)}>I^{(N)}$ is also evident.

4 The Kohn-Sham (KS) DFT

HK theorem allows the total energy $E[\rho_{0}]$ to exactly determine the ground state density $\rho_{0}$ via the variational Eq. [30], which is sufficient to determine the ground state properties. But, no direct expression of $E[\rho_{0}]$ \textcolorblackcan be found using HK theorem. In KS scheme, local single particle effective potential $v_{eff}({\bf r})=v_{ext}({\bf r})+v_{ee}({\bf r})$ is found which corresponds to the exact ground state density $\rho_{0}({\bf r})$ of the interacting system. This requires an assumption that the $v$ -representable ground state density of interacting system $\rho_{0}({\bf r})$ is also $v$ -representable in a non-interacting system $S$ . Once single particle effective potential is found the coupled equation for $N$ interacting electrons transforms to a single particle equation, called the KS equation. \textcolorblackThe eigenstates found as the solution of KS equation are the KS spin-orbitals. In ground state (GS), the $N$ electrons occupy the lowest lying $N$ spin-orbitals. Remember that, each energy eigenstate is doubly degenerate for up-down spin states. So, from total energy, using KS scheme, we can find a single particle Schrödinger equation as:

	$\displaystyle[-\frac{1}{2}\nabla^{2}+v_{eff}({\bf r})]\psi_{i}([\rho_{0}],{\bf x})=\epsilon_{i}\psi_{i}({\bf x})~{}~{}~{}\text{where, }~{}~{}~{}[\epsilon_{1}\leqslant\epsilon_{2}\leqslant...]$		(42)
	$\displaystyle\text{The kinetic energy:}~{}T^{S}[\rho_{0}]=\sum_{i=1}^{N}\int\psi^{\ast}_{i}([\rho_{0}],{\bf x})(-\frac{1}{2}\nabla^{2})\psi_{i}([\rho_{0}],{\bf x})d{\bf x}$		(43)
	$\displaystyle\text{The GS density:}~{}\rho_{0}({\bf r})=\sum_{i=1}^{N}\int\psi^{\ast}_{i}([\rho_{0}],{\bf x})\psi_{i}([\rho_{0}],{\bf x})d{\bf x}$		(44)

As, $T^{S}[\rho]$ is the kinetic energy of $S$ system, so, in the KS e-e interaction energy $E^{KS}_{ee}[\rho_{0}]$ , correlation-kinetic effect also contributes along with the combined effect (already present in xc) of Pauli principle and Coulomb repulsion. However, KS-DFT does not describe their individual contribution in $E^{KS}_{ee}[\rho_{0}]$ , whereas, in Q-DFT framework these terms are exactly defined (see, article 3.4 in (sahni2016, )). Since, the Coulomb part of $e-e$ interaction energy is the Hartree energy ( $E_{H}$ ) of Eq. [29], we can write $E^{KS}_{ee}[\rho_{0}]$ as:

	$\displaystyle E^{KS}_{ee}[\rho_{0}]=E_{H}[\rho_{0}]+E_{xc}^{KS}[\rho_{0}]~{}~{}\&~{}~{}v_{ee}({\bf r})=v_{H}({\bf r})+v_{xc}({\bf r})$		(45)
	$\displaystyle v_{H}([\rho_{0}],{\bf r})=\frac{\partial E_{H}[\rho_{0}]}{\partial\rho_{0}({\bf r})}~{};~{}~{}v_{xc}({\bf r})=v_{x}({\bf r})+v_{c}({\bf r})=\frac{\partial E_{xc}^{KS}[\rho_{0}]}{\partial\rho_{0}({\bf r})}$		(46)
	$\displaystyle v_{eff}=v_{ext}+v_{ee}=v_{ext}+v_{H}+v_{xc}=v_{ext}+v_{H}+v_{x}+v_{c}$		(47)

Finally, the Kohn-Sham total energy using the definitions in Eq. [43] and Eq. [29] :

\displaystyle E^{KS}=T^{S}[\rho_{0}]+\int\rho_{0}({\bf r})v_{ext}({\bf r})d{\bf r}+E_{H}[\rho_{0}]+E_{xc}^{KS}[\rho_{0}]

(48)

This $E_{xc}^{KS}[\rho_{0}]$ is not explicitly known and requires further approximations as described in Section 5.

4.1 Generalised KS Schemes

We have already discussed that the $\Psi([\rho_{0}],{\bf x})$ can be found from $\rho_{0}({\bf r})$ through constrained search formalism, as, there is no unique $B^{-1}$ mapping available (Please look into Sec. 3). This can be done by minimizing the universal density functional $F^{S}$ (of Eq. [37]) using Slater type trial wavefunction $\Psi^{S,t}$ that corresponds to $\rho_{0}({\bf r})$ . There can be infinite number of such trial wavefunctions. In generalised KS scheme, one searches for the infimum value of the expectation $\langle(\hat{T}+\hat{V}_{ee})\rangle$ with respect to all $N$ dimensional single Slater wavefunctions $\Psi^{S,t}$ that yields the ground state density $\rho_{0}({\bf r})$ .

F^{S}[\rho_{0}]=\inf_{\Psi^{S,t}\rightarrow\rho_{0}({\bf r})}\langle\Psi^{S,t}|(\hat{T}+\hat{V}_{ee})|\Psi^{S,t}\rangle

(49)

As, $\hat{V}_{ee}$ is not explicitly known in $S$ system, in standard KS scheme $\langle\hat{T}\rangle$ is used instead of $\langle(\hat{T}+\hat{V}_{ee})\rangle$ . Also, this allows to use single Slater wavefunction which is for non-interacting systems. Other two examples of generalised KS (GKS) schemes are Hartree-Fock KS scheme and screened non-local exchange (sx-LDA) scheme, those take a part of $\hat{V}_{ee}$ into account while minimizing (seidl1996, ):

		$\displaystyle\text{KS: }~{}~{}~{}~{}~{}~{}F^{S}[\rho_{0}]=\inf_{\Psi^{t}\rightarrow\rho_{0}({\bf r})}\langle\Psi^{t}\|\hat{T}\|\Psi^{t}\rangle$		(50)
		$\displaystyle\text{HF KS: }~{}~{}~{}~{}~{}~{}F^{S}[\rho_{0}]=\inf_{\Psi^{t}\rightarrow\rho_{0}({\bf r})}\left(\langle\Psi^{t}\|\hat{T}\|\Psi^{t}\rangle+E_{H}[\{\psi_{i}\}]+E_{x}[\{\psi_{i}\}]\right)$
		$\displaystyle\text{sx-LDA: }~{}~{}~{}~{}F^{S}[\rho_{0}]=\inf_{\Psi^{t}\rightarrow\rho_{0}({\bf r})}\left(\langle\Psi^{t}\|\hat{T}\|\Psi^{t}\rangle+E_{H}[\{\psi_{i}\}]+E^{sx}_{x}[\{\psi_{i}\}]\right)$

Here, the exchange energy $E_{x}[\{\psi_{i}\}]$ is not the general HF exchange and $E^{sx}_{x}[\{\psi_{i}\}]$ is a screened exchange interaction.

4.2 KS DFT for Fractional $N$ ; Derivative Discontinuity and Bandgap:

With the above background, we now move on to discuss the origin of the underestimation of bandgap in KS scheme. The Euler-Lagrange equations (Eq. [38]) for KS equivalent $S$ system having fractional $N\pm f$ number of electrons are expressed using Eq. [46-48] as (perdew1997, ):

$\displaystyle\left.\frac{\delta T^{S}}{\delta\rho}\right\|_{N+f}+v_{ext}+v_{H}+v^{N+f}_{xc}$	$\displaystyle=-A^{(N)}$
$\displaystyle(-)\phantom{12}\left.\frac{\delta T^{S}}{\delta\rho}\right\|_{N-f}+v_{ext}+v_{H}+v^{N-f}_{xc}$	$\displaystyle=-I^{(N)}$
$\displaystyle\left.\frac{\delta T^{S}}{\delta\rho}\right\|_{N+f}-\left.\frac{\delta T^{S}}{\delta\rho}\right\|_{N-f}+v^{N+f}_{xc}-v^{N-f}_{xc}$	$\displaystyle=I^{(N)}-A^{(N)}$	(51)

It is proved for $N=1$ and argued for $N>1$ (sagvolden2008, ) that the KS energy bandgap ( $E_{g}^{KS}$ ), i.e., the HOMO-LUMO (or, VBM-CBM) gap of the $S$ system arises due to the discontinuity of derivative of kinetic energy (sham1983, ), i.e.

$\displaystyle E_{g}^{KS}=\left.\frac{\delta T^{S}}{\delta\rho}\right\|_{N+f}-\left.\frac{\delta T^{S}}{\delta\rho}\right\|_{N-f}$	$\displaystyle=\epsilon^{(N)}_{N+1}-\epsilon^{(N)}_{N}=\epsilon^{(N)}_{LUMO}-\epsilon^{(N)}_{HOMO}$	(52)
$\displaystyle\therefore I^{(N)}-A^{(N)}=\epsilon^{(N)}_{N+1}-\epsilon^{(N)}_{N}$	$\displaystyle+v^{N+f}_{xc}-v^{N-f}_{xc}\Rightarrow E_{g}=E_{g}^{KS}+\Delta v_{xc}$
$\displaystyle\&~{}~{}~{}\Delta v_{xc}=v^{N+f}_{xc}-v^{N-f}_{xc}$	$\displaystyle=\left(I^{(N)}-A^{(N)}\right)-\left(\epsilon^{(N)}_{N+1}-\epsilon^{(N)}_{N}\right)$	(53)

For fractional electron number $(N-f)$ , the spin-orbitals $\psi_{i}({\bf x}_{i})$ of Eq. [42] are occupied completely till $(N-1)$ and the $N$ -th spin-orbital should have fractional occupancy of fraction $f$ . The ground state energy density $\rho_{0}({\bf r})$ , and, the total energy in KS equivalent $S$ system for $(N-f)$ electrons follow the ensemble equation for fractional electron (Eq. [39]), so that:

	$\displaystyle\rho^{(N-f)}=\sum_{i}^{N-1}\|\psi_{i}([\rho_{0}],{\bf x}_{i})\|^{2}+f\|\psi_{N}([\rho_{0}],{\bf x}_{N})\|^{2}=(1-f)\rho_{0}^{(N-1)}+f\rho_{0}^{(N)}$
	$\displaystyle E^{(N-f)}=(1-f)E^{(N-1)}+fE^{(N)}\Rightarrow\left.\frac{\partial E}{\partial N}\right\|_{N-}=E^{(N)}-E^{(N-1)}=\epsilon^{(N)}_{N}$		(54)

\textcolor

blackFrom Eq. [34], we see that $E^{(N)}-E^{(N-1)}=-I^{(N)}$ . So, we finally reach to the Ionization Potential (IP) theorem (perdew1997, ):

The highest occupied Kohn-Sham eigenvalue for $N$ electron system (even fractional) is the minus of the ionization energy of the system ; $\epsilon^{(N)}_{N}=-I^{(N)}$ .

5 Jacob’s Ladder of Exchange-Correlation Functional

The xc potential $v_{xc}$ or energy $E_{xc}$ \textcolorblackcan not be exactly calculated for many-electron system. There are several levels of approximations used to express the exchange correlation of KS system. It is like climbing a ladder; more rungs one climbs, more accuracy is achieved. Obviously, to achieve that, more effort is needed, in terms of mathematical complexcity, and, computational costing. From now if spin-orbitals are referred as $\psi({\bf x})$ or orbitals as $\phi({\bf r})$ they will represent the density dependent versions of KS scheme $\psi([\rho_{0}],{\bf x})$ and $\phi([\rho_{0}],{\bf r})$ , if not explicitly defined otherwise.

[scale=0.5]Jacob_s-ladder

Figure 2: Jacob’s Ladder in DFT.

The first rung of Jacob’s ladder (Fig. 2) starts with the introduction of the easiest model for calculating xc, the homogeneous electron gas (HEG) model or Jellium model (perdew2001, ). In this model, the non-relativistic electron cloud is taken to be homogenous, i.e., density is everywhere same $\left(\rho_{0}({\bf r})=const.\right)$ . A solid state system can be assumed to be a box (Wigner cell) with periodic boundary condition and the electron orbitals can be described by plane waves $\phi_{i}=\frac{1}{\sqrt{\mathscr{V}_{cell}}}e^{\sqrt{-1}{\bf k}.{\bf r}}$ , where $\mathscr{V}_{cell}$ is the volume of the cell. As, $k=\frac{2\pi}{\lambda}$ and electron wavefunction should vanish at cell boundary $\pm L/2$ , for one dimension, we can understand that one electron should occupy $2\pi/L$ length in real space. Equivalently, in three dimension it should occupy a volume $\frac{(2\pi)^{3}}{\mathscr{V}_{cell}}$ . Hence, number of electrons in the wave vector range ${\bf k}$ to ${\bf k}+d{\bf k}$ , considering both spins, is given by $2\frac{\mathscr{V}_{cell}}{(2\pi)^{3}}d{\bf k}$ and $d{\bf k}=4\pi k^{2}dk$ . The total number of electron $N$ , and, the electron density $\rho$ can be found by integrating $dk$ till the Fermi wave vector $k_{F}$ corresponding to the highest occupied energy, i.e., Fermi energy $E_{F}$ .

\displaystyle N=\int_{0}^{k_{F}}2\frac{\mathscr{V}_{cell}}{(2\pi)^{3}}4\pi k^{2}dk=\frac{\mathscr{V}_{cell}k_{F}^{3}}{3\pi^{2}}~{};~{}~{}\rho_{0}=\frac{k_{F}^{3}}{3\pi^{2}}\Rightarrow k_{F}=\{3\pi^{2}\rho_{0}\}^{1/3}\propto\rho_{0}^{1/3}

(55)

5.1 Local Density Approximation (LDA)

In local density approximation, the electron gas is considered to be locally uniform ( $\rho_{0}({\bf r})=\rho_{0}$ ) as in HEG, and, the exchange energy is expressed in term of exchange energy per particle $\varepsilon_{x}^{HEG}$ or in term of exchange energy density $\tilde{\varepsilon}_{x}^{HEG}=\varepsilon_{x}^{HEG}\rho_{0}$ (for derivation of $\varepsilon_{x}^{HEG}$ value, see Sec. 4.3 of Ref. (engel2013, )) as:

	$\displaystyle E_{x}^{LDA}[\rho]=\int\varepsilon_{x}^{HEG}[\rho_{0}]~{}\rho_{0}d{\bf r}~{};~{}\varepsilon_{x}^{HEG}[\rho_{0}]=-\frac{3k_{F}({\bf r})}{4\pi}=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho_{0}^{1/3}$		(56)
	$\displaystyle\Rightarrow E_{x}^{LDA}[\rho]=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\int\rho_{0}^{4/3}d{\bf r}~{}~{}\&~{}~{}v_{x}^{LDA}[\rho_{0}]=-\left(\frac{3}{\pi}\right)^{1/3}\rho_{0}^{1/3}=-\frac{k_{F}}{\pi}$		(57)

If we take a sphere where one electron is present, then the density can be presented as $\rho_{0}=\frac{1}{4/3\pi r_{w}^{3}}$ . Such sphere is called Wigner-Seitz sphere of radius $r_{w}$ .

r_{w}=\left(\frac{3}{4\pi\rho_{0}}\right)^{1/3}=\left(\frac{9\pi}{4}\right)^{1/3}\frac{1}{k_{F}}~{}~{}\&~{}~{}~{}\varepsilon_{x}^{HEG}[\rho_{0}]=-\frac{3}{4\pi}\left(\frac{9\pi}{4}\right)^{1/3}\frac{1}{r_{w}}

(58)

The correlation is exactly known for two extreme limits, high density weak coupling limit ( $r_{w}\rightarrow 0$ ) and low density strong coupling limit ( $r_{w}\rightarrow\infty$ ).

	$\displaystyle\left.\varepsilon_{c}\right\|_{r_{w}\rightarrow 0}$	$\displaystyle=$	$\displaystyle c_{0}\ln r_{w}-c_{1}+c_{2}r_{w}\ln r_{w}-c_{3}r_{w}+...$
	$\displaystyle\left.\varepsilon_{c}\right\|_{r_{w}\rightarrow\infty}$	$\displaystyle=$	$\displaystyle-\frac{d_{0}}{r_{w}}+\frac{d_{1}}{r_{w}^{3/2}}-\frac{d_{2}}{r_{w}^{2}}+...$		(59)

For high density limit, the leading contribution can be calculated using RPA, and, the next higher level contribution is calculated using first order response function (second order exchange). The first two positive constants are $c_{0}=0.031091$ and $c_{1}=0.0467$ . A complete RPA calculation is provided by von Barth and Hedin, as well as, by Vosko, Wilk and Nusair (von1972, ; vosko1980, ). For, low density limit, the first two constants $d_{0}$ and $d_{1}$ can be estimated from the (Madelung) electrostatic energy and the zero point vibrational energies of the Wigner crystal. For real system having intermediate range of densities, the simplest approach is to interpolate in-between, and, there are several methods developed over years (vosko1980, ; perdew1981, ). A systematic discussion on the development of LDA exchange is presented in Chap-6 of (dreizler2012, ). Generalization of LDA for spin polarized system is known as LSDA, and, the exchange is $\varepsilon^{HEG}[\rho_{0\uparrow},\rho_{0\downarrow}]$ . As the HEG model better suits solids than molecules, LSDA performs better in solids. The bond lengths and lattice constants \textcolorblackpredicted by LSDA are manageable. But in the case of band gap or ionization energy estimation, the scenario is far worse.

To understand the reason, let us go back to the exchange correlation potential in Eq. [57]. As soon as the density $\rho_{0}({\bf r})$ starts to fall exponentially as $\rho_{0}({\bf r})\sim e^{-\alpha r}$ , the exchange potential $v_{x}^{LDA}[\rho_{0}]\sim e^{-\alpha r/3}$ . In asymptotic region ( $r\rightarrow\infty$ ) the the exponential fall of exchange potential violates the condition $\lim\limits_{r\rightarrow\infty}v_{x}\sim\frac{1}{r}$ (levy1984, ). Similar problem is experienced for correlation which should follow the power law nature (almbladh1985, ). Another factor of the failure of LDA is termed as Delocalization error which is discussed in [B3].

5.2 Generalised Gradiant Approximation (GGA)

LDA uses homogeneous electron gas model, but electron density can never be homogenous in any real system. As a starting point HEG model is good, but then, one have to introduce next level of corrections. The general choice is the slowly varying HEG model for which $\varepsilon_{x}$ can be calculated using gradient expansion approximation (GEA) for high density limit as:

E_{x}^{GEA}=\int\varepsilon_{x}^{HEG}[\rho_{0}]~{}\rho_{0}\left[1+\tilde{\mu}s^{2}+\cdots\right]d{\bf r}~{};~{}~{}s=\frac{|\nabla\rho_{0}|}{2k_{F}\rho_{0}}=\frac{|\nabla\rho_{0}|}{2\{3\pi^{2}\rho_{0}\}^{1/3}\rho_{0}}

(60)

Here $\{\tilde{\mu},\cdots\}$ are constants which can be calculated from GEA. Antoniewicz and Kleinman have calculated $\tilde{\mu}^{GEA}=10/81$ (antoniewicz1985, ). Though, the mathematical derivation is rigorous for GEA, however, for atoms and molecules, the slowly varying idea can not be justified. But, it has given a clue, and, finally Perdew and Wang (PW91) (PW91, ) have made the generalised gradient approximation simple by introducing analytical parameters. \textcolorblackFurther improvements to this have lead to the Perdew, Bruke, Ernzerhof (PBE) form (PBE, ) of GGA. Using exchange enhancement factor $\tilde{f}(s)$ , $E_{x}^{GGA}$ becomes:

	$\displaystyle E_{x}^{GGA}=\int\varepsilon_{x}^{HEG}[\rho_{0}]~{}\rho_{0}\tilde{f}^{\prime}(s)d{\bf r}=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\int\rho_{0}^{4/3}\tilde{f}(s)d{\bf r}$		(61)
	$\displaystyle\tilde{f}^{PBE}(s)=1+\frac{\tilde{\mu}s^{2}}{1+\tilde{\mu}s^{2}/\kappa}~{};~{}~{}\tilde{\mu}=0.21951\approx 1.78\tilde{\mu}^{GEA},\kappa\leq 0.804$		(62)

The exchange enhancement factor, which is generally expressed as a power series of $s$ , recovers the uniform gas limit for $\tilde{f}(s)=1$ . The $\tilde{\mu}\approx 2\tilde{\mu}^{GEA}$ is followed by PW91, PBE and other GGA functionals and is appropriate for expressing exchange energies of neutral atoms. But in solids the densities are almost slowly varying, so $\tilde{\mu}=\tilde{\mu}^{GEA}$ is a better choice. Thus, in the modified version of PBE for solids (PBEsol) $\tilde{\mu}$ is taken as $\mu^{GEA}$ (PBEsol, ).

The density falls exponentially as $\rho_{0}({\bf r})\sim e^{-\alpha r}$ ( $\alpha$ is decay constant) and $\alpha$ can be expressed in term of chemical potential $\mu_{-}^{(N)}$ (do not get confused with $\tilde{\mu}$ of Eq. [60]) (vLB, ; levy1984, ). Thus, Eq. [40], $\mu_{-}^{(N)}=-I^{(N)}$ , and, the IP theorem $\epsilon^{(N)}_{N}=-I^{(N)}$ yield:

\rho_{0}({\bf r})=Ne^{-\alpha r},\text{ where, }\alpha=2\sqrt{-2\mu_{-}^{(N)}}=2\sqrt{2I^{(N)}}=2\sqrt{-2\epsilon^{(N)}_{N}}

(63)

So, density at asymptotic region depends on the HOMO eigenvalue $\epsilon^{(N)}_{N}$ . Now, $\rho_{0}$ from Eq. [63] is inserted into the Eq. [61] to finds $E_{x}^{GGA}$ , and, the corresponding potential $v_{x}^{GGA}=\frac{\partial E_{x}^{GGA}}{\partial\rho_{0}({\bf r})}$ . \textcolorblackThus, $\tilde{f}^{\prime}(s)$ can be found.

The $v_{x}^{GGA}$ should follow the asymptotic limit $\lim\limits_{r\rightarrow\infty}v_{x}\sim-\frac{1}{r}$ . This is the starting idea of the semi-empirical GGAs. Furthermore, Savin et al. have showed that for finite systems, the accurate treatment of exchange-correlation potential, particularly in the asymptotic region, leads to eigenvalue differences between the HOMO and LUMO energies close to true bandgap (savin1998, ). This is the philosophy of semi-empirical functional which has been first introduced by Becke (Becke86) (becke1986, ). Parametrization of Becke86 and PW91 potentials have also been proposed for producing better asymptotic behaviour (harbola2002, ).

Similar philosophy of producing exact nature of exchange potential at asymptotic region for better electronic structure prediction has been followed by vanLeeuwen and Baerends (vLB) and the corrected LDA exchange-correlation potential (also known as LB94) has been built as (vLB, ):

v_{xc}^{vLB}({\bf r})=[v_{x}^{LDA}({\bf r})+v_{x}^{vLB}({\bf r})]+v_{c}^{LDA}({\bf r})~{};~{}~{}v_{x}^{vLB}({\bf r})=-\frac{\beta z^{2}\rho^{1/3}({\bf r})}{1+3\beta\ z\sinh^{-1}(z)}

(64)

Here, $z=\frac{|{\nabla}\rho({\bf r})|}{\rho^{4/3}({\bf r})}=\frac{s}{2\{3\pi^{2}\}^{1/3}}$ and $\beta=0.05$ , It has been employed to study the effect of the correct asymptotic behaviour of potential on the response properties of atoms (banerjee1999, ), and, is a very successful potential in describing different properties (e.g, photoionization, absorption spectra, etc.) in atomic and molecular physics due to its correct asymptotic behaviour (stener2000, ; pertot2017, ). The vLB corrected LDA has been first applied to solids by Singh et al. (singh2013, ; singh2015, ). In its original form, the $\beta$ has been taken as constant for every material, which may not be always logical. In (singh2016, ), Singh et al. have varied the $\beta$ using the IP theorem ( $\epsilon^{(N)}_{N}=-I^{(N)}$ ), and, have applied to diverse set of materials, which indeed becomes a successful method for proper bandgap prediction. In this IP-vLB work, for optimal choice, $\beta$ has been varied until the highest-occupied molecular orbital (HOMO) eigenvalue $\epsilon_{N}^{(N)}$ matches with the minus of ionization energy ( $-I^{(N)}$ ).

The material dependency can also be availed through the muffin tin radius in localised orbital schemes, or, in better form in full potential NMTO (FP-NMTO) method through a combination of hard sphere and muffin tin radii (NMTO, ; NMTO1, ). The minimal basis set used in FP-NMTO method is more accurate and flexible than LMTO-ASA. Being a full-potential method, it can handle the interstitial region more accurately. The larger overlap between muffin tin orbitals in FP-NMTO can produce the curvature of potential in the region in-between atoms can be produced more precisely. Also, the introduction of higher order energy derivatives in the process of removing the energy dependency of the basis set minimizes the error. This has been utilised by Datta et al. by incorporating the original vLB potential in self-consistent FP-NMTO, which produces almost exact bandgaps for traditional group IV (Si, Ge) and group III-V (GaAs, InP, etc.) semiconductors (datta2019, ). This method is as fast as any other GGA calculation. So, this type of potential-only correction which is not coming from any functional derivative of energy can also provide better electronic structure prediction with tricky tuning as in (singh2016, ; datta2019, ).

A comparison of LDA, PBE and vLB xc potential in $[100]$ plane for Ge, calculated using FP-NMTO method, is presented in Fig. 3 (datta2019, ). The LDA potential is smooth, whereas, PBE is showing a variation in the interstitial region (region in-between Ge atoms) and vLB produces larger variation. These variations are is due to the involvement of the gradient correction of charge density.

Figure 3: FP-NMTO XC potential (z-axis is range) for LDA, PBE and LDA+vLB in

\left[100\right]

-plane.

Direct potential correction has been followed by other methods, as in Becke Johnson (BJ) (becke2006, ) GGA. Being a derivative of energy functional, the potential can only be defined to within an arbitrary constant, so, $\lim\limits_{r\rightarrow\infty}v_{x}\sim-\frac{1}{r}+const.$ is also possible. BJ have proposed an effective exchange potential, which is not derived from energy derivative, followed that. They have shifted the exchange potential so that HOMO energy equals to its exact Hartree-Fock HOMO energy.

Armiento and Kümmel (AK13) have proposed an exchange enhancement factor $\tilde{f}^{AK13}$ , so that, the potential can be found by functional derivative, and then, like the BJ scheme, the asymptotic constant shift of exchange potential has been utilised (armiento2013, ).

\tilde{f}^{AK13}(s)=1+B_{1}s\ln(1+s)+B_{2}s\ln(1+\ln(1+s));~{}v_{x}^{mod}({\bf r})=v_{x}({\bf r})-\lim\limits_{r\rightarrow\infty}v_{x}({\bf r})

(65)

Where $B_{1}=3/5\tilde{\mu}^{GEA}+8\pi/15$ and $B_{2}=\tilde{\mu}^{GEA}-B_{1}$ . In this process, the HOMO eigenvalue is also shifted as $\alpha=2\sqrt{-2\left(\epsilon^{(N)}_{N}-\lim\limits_{r\rightarrow\infty}v_{x}({\bf r})\right)}$ in Eq. [63].

Later, they have pointed out some drawbacks of such nonvanishing asymptotic exchange potentials, those demand constant shift, though, the proper production of derivative discontinuity made these successful in reproducing atomic and molecular properties, and, atomic-shell structures. Bandgaps predicted by these are close to experimental values establishing their usefulness in semiconductor theory (aschebrock2017, ).

Gritsenko et al. (GLLB) have noticed from the optimised potential method (OPM) result that, the shell closure like nature of exchange potential can be reproduced by a combination of two terms; one is a screening type Slater exchange (Eq. [26]]) originating from Fermi exchange hole, and, another is the response of exchange only pair-correlation function $g_{x}({\bf r},{\bf r}^{\prime})$ of Eq. [26], written as (gritsenko1995, ):

	$\displaystyle v_{x}({\bf r})$	$\displaystyle=v_{x}^{S}({\bf r})+v_{x}^{resp}({\bf r})$		(66)
	$\displaystyle\text{where,}~{}~{}v_{x}^{resp}({\bf r})$	$\displaystyle=\frac{1}{2}\int\rho_{0}({\bf r}^{\prime})\left[\int\frac{\rho_{0}({\bf r}^{\prime\prime})}{\|{\bf r}^{\prime}-{\bf r}^{\prime\prime}\|}\frac{\delta g_{x}({\bf r}^{\prime},{\bf r}^{\prime\prime})}{\delta\rho_{0}({\bf r})}d{\bf r}^{\prime\prime}\right]d{\bf r}^{\prime}$		(67)

GLLB have approximated $v_{x}^{resp}({\bf r})$ using HOMO energy $\epsilon^{(N)}_{N}$ and eigenstates $\{\psi_{i}\}$ as:

v_{x}^{resp,GLLB}({\bf r})\simeq K_{g}\sum_{i}^{N}\sqrt{\epsilon^{(N)}_{N}-\epsilon^{(N)}_{i}}\frac{|\psi_{i}({\bf x})|^{2}}{\rho_{0}({\bf r})}

(68)

where, $K_{g}$ is calculated from HEG model as $8\sqrt{2}/3\pi^{2}$ . Kuisma et al. have combined PBEsol correlation with this GLLB exchange to produce GLLB-SC $xc$ potential aimed for solids (kuisma2010, ):

\displaystyle v_{xc}^{GLLB-SC}({\bf r})=2\varepsilon_{x}^{PBEsol}({\bf r})+K_{g}\sum_{i}^{N}\sqrt{\epsilon^{(N)}_{N}-\epsilon^{(N)}_{i}}\frac{|\psi_{i}({\bf x})|^{2}}{\rho_{0}({\bf r})}+v_{c}^{PBEsol}({\bf r})

(69)

These GLLB and GLLB-SC potentials produce the derivative discontinuity for solids, because, a small $\delta\rho$ addition to density corresponding integer electron opens up a new orbital with fractional occupation and $\epsilon^{(N)}_{N}$ jumps. GLLB-SC performs better than GLLB due to incorporation of PBEsol GGA correlation (tran2018, ).

{svgraybox}

B3: Delocalization Error

Reproduction of piecewise linearity of total energy of HK systems (see [B2]) is a necessary condition to follow by every (semi-)local approximations. For system with fractional charge, LDA energy ( $E$ ) curve with electron number ( $N$ ) gives a piecewise convex nature. It means these functionals should calculate very low energy for fractional charge. As fractional charge arises due to electron exchange between atoms, it signifies the electron delocalization, and, functionals with convex nature wrongly delocalizes electron. So, this error is termed as Delocalization error (mori2008, ). This piecewise convex behaviour is present in many semi-local functionals. On the other hand, the HF energy curve is piecewise concave, and, imposing localization of xc kernel through self-interaction correction (SIC), e.g., Perdew-Zunger PZ-SIC, follows similar wrong concave nature (see, left panel of Fig. 5) (mori2008, ; vydrov2007, ). Long range corrected hybrid functional LC- $\omega$ PBE produces almost exact piecewise linear character. Besides these numerical evidences, recently the concavity of HF energy is proved mathematically (li2017, ). However, the piecewise convexity of semi-local functionals is not guaranteed for all fractions, especially for low fractions, as seen in the right panel of Fig. 5.

5.3 Meta-GGA

LDA starts with the homogeneous distribution of density and GGA introduces the \textcolorblackeffect of first order spatial variation of density ( $\nabla\rho_{0}$ ) \textcolorblackwithin xc potential. The next level of correction should naturally come through the introduction of the second order term $\nabla^{2}\rho_{0}$ . Some functions have been proposed, but, there is a problem with such terms; the corresponding potential becomes fourth order gradient of density, which blows out near the nuclei making the calculation highly unstable. The kinetic energy density term $\tau({\bf r})$ provides an alternative pathway.

\tau({\bf r})=\sum_{\sigma}\tau_{\sigma}({\bf r})=\frac{1}{2}\sum_{\sigma}\sum_{i}^{occ}|\nabla\psi_{i}({\bf r},\sigma)|^{2}

(70)

The gradient expansion of kinetic energy provides a flexibility, we can add $\nabla^{2}\rho_{0}$ dependent term in kinetic energy density not changing the kinetic energy itself. This is because $\int\nabla^{2}\rho_{0}d{\bf r}$ vanishes (Gauss’ theorem) if integrated over any finite system. Using gradient expansion (may refer to Sec. 4.4.3 of (engel2013, )) of kinetic energy density upto second order, one can reach to the expansion (brack1976, ):

\tau^{GEA}({\bf r})=\tau^{HEG}+\frac{1}{72}\frac{|\nabla\rho_{0}({\bf r})|^{2}}{\rho_{0}({\bf r})}+\frac{1}{6}\nabla^{2}\rho_{0}({\bf r})~{}~{};~{}~{}\tau^{HEG}=\frac{3(6\pi^{2})^{2/3}}{10}\rho_{0}^{5/3}({\bf r})

(71)

So, kinetic energy density $\tau({\bf r})$ can serve as a substitute of $\nabla^{2}\rho_{0}({\bf r})$ , and, the story of meta-GGA (MGGA) begins. In MGGA, the exchange enhancement factor is defined as: $\tilde{f}_{x}^{MGGA}\equiv\tilde{f}_{x}(\rho_{0},\nabla\rho_{0},\tau)$ .

Perdew, Kurth, Zupan and Blaha (PKZB) have further used the idea of PBE GGA to develop a semi-empirical meta-GGA satisfying spin-scaling, uniform density scaling and Lieb-Oxford lower bound constraints (PZKB, ). \textcolorblackA list of constraints used in DFT approximations can be found in Ref. (SCAN, ). In reduced functional of $p$ and $\tilde{q}$ , $\tilde{f}_{x}$ is expressed as:

	$\displaystyle\tilde{f}^{PKZB}_{x}(p,\tilde{q})=1+\frac{\zeta}{1+\zeta/\kappa}~{}~{}\text{where, }~{}~{}p=s^{2},~{}\tilde{q}$	$\displaystyle=\frac{3\tau}{2(3\pi^{2})^{2/3}\rho_{0}^{5/3}}-\frac{9}{20}-\frac{p}{12}$		(72)
	$\displaystyle\zeta^{PZKB}=\frac{10}{81}p+\frac{146}{2025}\tilde{q}^{2}-\frac{73}{405}p\tilde{q}$	$\displaystyle+\left[D+\frac{1}{\kappa}\left(\frac{10}{81}\right)^{2}\right]p^{2}$		(72)

The first two term in $\zeta^{PZKB}$ do not contain $\kappa$ , so that the exchange enhancement factor recovers the slowly varying limit of gradient expansion (GEA) of enhancement factor up to the fourth order of $\nabla\rho_{0}({\bf r})$ as:

\tilde{f}_{x}^{GEA}=1+\frac{10}{81}p+\frac{146}{2025}\tilde{q}^{2}-\frac{73}{405}p\tilde{q}+Dp^{2}+\cdots

(73)

In PKZB scheme, the empirical parameter $D=0.113$ is estimated by minimizing the mean absolute error of the atomization energies for a set of molecules. PKZB functional successfully determines the surface and atomization energies but it overestimates the bond lengths. Also, it remains unsuccessful in describing hydrogen bonded systems.

Tao et al. have attributed this failure in bond length prediction to the PKZB exchange, and, have proposed a corrected exchange called TPSS (Tao-Perdew-Staroverov-Scuseria) (TPSS, ). TPSS is the first non-empirical MGGA which makes the empirical parameter of PKZB, $D=0$ . Further it includes an additional constraint: making exchange potential finite at the nucleus for the ground state of one and two electronic system. The single orbital kinetic energy density reduces to the Weizsäcker functional (for spin polarised system, single orbital means two spin channels) $\tau^{W}$ and two dimensionless quantities are useful (using $\tau^{HEG}$ of Eq. [71]):

\tilde{z}=\frac{\tau^{W}}{\tau}~{};~{}~{}\tilde{\alpha}=\frac{\tau-\tau^{W}}{\tau^{HEG}}~{};~{}~{}\text{where, }\tau^{W}=\frac{|\nabla\rho_{0}({\bf r})|^{2}}{8\rho_{0}({\bf r})}

(74)

The iso-orbital indicator $\tilde{\alpha}$ differentiate between different local bonding environments: $\tilde{\alpha}=0$ or $\tilde{z}=1$ represents the covalent single bond, $\tilde{\alpha}\approx 1$ i.e., $\tau-\tau^{W}\approx\tau^{HEG}$ indicates that after bonding there are enough electron cloud available to replicate HEG, so, represents metallic bonding, and, $\tilde{\alpha}\gg 1$ represents weak bonding (sun2013, ). For two electron system ( $\tilde{\alpha}=0$ ), the PKZB-MGGA enhancement factor of Eq. [72] reduces to GGA form, like Eq. [62]. The corresponding exchange potential, which has a $\nabla^{2}\rho_{0}({\bf r})$ term, diverges at nucleus. To avoid that, the constraint $\left.\frac{d\tilde{f}_{x}}{ds}\right|_{s=0.376}=0$ has been imposed in TPSS scheme, where, $s=0.376$ is the value of $s$ at nucleus for two-electronic system. Satisfying all these constraints, they have reconstructed the $\zeta$ of Eq. [72] as a complicated analytical functional: $\zeta[p,\tilde{q},\tilde{z}]$ (TPSS, ).

TPSS functional has improved the \textcolorblackefficiency of the estimation of lattice constants for solids as well as the bond lengths for molecules, hydrogen bonded complexes whereas, the mean error in calculation of bulk-moduli of solids increases with respect to PKZB functional. It is less accurate than the PBE GGA for calculation of critical pressure of structural phase transition of solids as well.

To solve such problems, Sun et al. have introduced a functional which satisfies a larger set of constraints and appropriate norms. They have termed it as strongly constrained and appropriately normed (SCAN) functional (SCAN, ). In construction of SCAN functional, they have interpolated $\tilde{f}^{PKZB}$ in-between covalent bonding case ( $\tilde{\alpha}=0$ ), and, metallic bonding case ( $\tilde{\alpha}=1$ ). \textcolorblackThen, they have further extrapolated for hydrogen bonded systems ( $\tilde{\alpha}\rightarrow\infty$ ). There is a switching function defined, which switches between $\tilde{\alpha}=0$ and $\tilde{\alpha}=1$ in the expression of $\tilde{f}^{SCAN}$ . SCAN functional is successful in describing molecular energies, barrier heights of chemical reactions, as well as, in predicting lattice constants, mechanical stability of solids (sun2016, ; zhang2018, ). However, SCAN functional significantly overestimates the magnetization of elemental ferromagnetic materials (Fe, Co, and Ni). In a recent study, the reason is attributed to the insensitivity of the switching function to $\tilde{\alpha}$ for some particular range, and, oversensitivity in another range (mejia2019, ). They have proposed a deorbatalized version of SCAN, called as SCAN-L where $\tilde{\alpha}$ is a function of $[\rho,\nabla\rho,\nabla^{2}\rho]$ .

We have seen that (Eq. [24] using Eq. [14] & 15) the Fermi exchange hole $\rho_{x}({\bf r},{\bf r}^{\prime})$ can be expressed in term of first order density matrix. The Fermi hole is highly delocalized. A density matrix expansion (DME) under general coordinate transformation can make the Fermi hole localized, which reduces the difficulty in modelling the conventional Fermi hole. Tao and Mo (TM) have followed similar type of coordinate transformation as done in (van1998, ), and, formulate a MGGA exchange enhancement factor $\tilde{f}_{x}^{TM-DME}$ (TAO-MO, ). In the HEG limit, the exchange energy functional and the Fermi hole from the DME are exact, however, this fails for slowly varying densities. To solve this, they have interpolated between the DME and slowly varying densities, as well as, between the corresponding enhancement factors. The slowly varying enhancement factor $\tilde{f}_{x}^{TM-SC}$ is found in accordance to the fourth order gradient expansion of $\tilde{f}_{x}^{GEA}$ (Eq. [73]). The TM interpolated enhancement factor using the $\tilde{z}$ dependent weight factor $w$ as:

\displaystyle\tilde{f}_{x}^{TM}=w~{}\tilde{f}_{x}^{TM-DME}+(1-w)\tilde{f}_{x}^{TM-SC}~{};~{}~{}w=\frac{\tilde{z}^{2}+3\tilde{z}^{3}}{(1+\tilde{z}^{3})^{2}}

(75)

So, TM have used $\tilde{z}$ as the iso-orbital indicator,which can differentiate between single-orbital ( $\tilde{z}=0$ ) and slowly varying $\tilde{z}\approx 1$ bondings. Another widely used indicator in MGGA is $\tilde{t}^{-1}=\frac{\tau}{\tau^{HEG}}$ , which can differentiate between covalent and non-covalent bonding, but, can not identify single-orbital regions. As mentioned earlier, TPSS and SCAN have used $\tilde{\alpha}$ , which is better than these two and can be directly related with electron localization function: $ELF=\frac{1}{1+\tilde{\alpha}^{2}}$ .

However, $\tilde{\alpha}$ dependent MGGA (e.g., SCAN) faces numerical instability originating from sharp oscillations in the xc potential, and, can only be eliminated with very fine grids (furness2019, ; bartok2019, ). To solve the problem, Furness and Sun have proposed a indicator $\tilde{\beta}=\frac{\tau-\tau^{W}}{\tau+\tau^{HEG}}$ (furness2019, ). Bartók and Yates have indicated regularised SCAN (rSCAN) functional, \textcolorblackwhere the switching function is expressed as $\tilde{\alpha^{\prime}}=\frac{\tilde{\alpha}^{3}}{\tilde{\alpha}^{2}+0.001}$ . They have used $\tau^{HEG}+10^{-4}$ instead of $\tau^{HEG}$ in Eq. [74] to regularize the $\tilde{\alpha}$ for only very small values (bartok2019, ). rSCAN improves numerical performance of SCAN at the expense of breaking constraints known from the exact xc functional. In SCAN, LSDA limit for HEG can be recovered through the constraint $\tilde{\alpha}\rightarrow 1$ , but, rSCAN loses the correct HEG description. The correct uniform- and non-uniform- scaling properties of $\alpha$ , as well as, the correct HEG limit of $E_{xc}$ are restored in r²SCAN (furness2020, ) using a regularization parameter $\eta$ as:

\tilde{\alpha}^{\prime\prime}=\frac{\tau-\tau^{W}}{\tau^{HEG}+\eta\tau^{W}}~{};~{}~{}\text{where, }\eta=10^{-3}

(76)

In MGGA regime, potential-only functionals (those are not functional derivative of exchange energy, like vLB in GGA section) are available as well. The Becke Johnson (BJ) semi-empirical potential is given as (becke2006, ):

v_{x,\sigma}^{BJ}=v_{x}^{BR}({\bf r})+\frac{1}{\pi}\sqrt{\frac{5}{6}}\frac{\tau({\bf r})}{\rho_{0}({\bf r})}

(77)

Here $v_{x}^{BR}$ is Becke Russel GGA potential (see, (becke1989, )). As the BJ potential have not been formed in the conventional way as any other MGGA functional, so, some veterans in this field do not want to consider it as MGGA tran2019 . Tran and Blaha have modified the BJ potential (TBmBJ) by introducing a density gradient dependent weight factor $w$ as (tran2009, ):

	$\displaystyle v_{x,\sigma}^{TBmBJ}=w~{}v_{x}^{BR}({\bf r})+$	$\displaystyle(3w-2)\frac{1}{\pi}\sqrt{\frac{5}{6}}\frac{\tau({\bf r})}{\rho_{0}({\bf r})}$		(78)
	$\displaystyle\text{where,}~{}~{}w=-0.012+1.023\sqrt{\mathscr{G}}~{}~{}$	$\displaystyle;~{}~{}\mathscr{G}=\frac{1}{{\mathscr{V}}_{cell}}\int_{cell}\frac{1}{2}\frac{\|\nabla\rho_{0}({\bf r}^{\prime})\|}{\rho_{0}({\bf r}^{\prime})}d{\bf r}^{\prime}$

Here, the constant $1.023$ is in (Bohr^1/2) unit. TBmBJ functional shows great improvement in bandgap prediction, and, optical properties calculations for semiconductors and insulators which will be discussed in the next section.

5.4 Hybrid Functional Method

The LDA, GGA , MGGA are the so called (semi-)local functionals, as the exchange energy density at any point ${\bf r}$ depends on the density, its first and second order gradients, and, atmost on the orbitals and/or their gradient, but, only at that point ${\bf r}$ . Whereas, for any real system, the exchange energy is expressed in term of density as in Eq. [27]. Due to the non-local nature of the exchange hole, the exchange energy integrand is fully non-local. This is why the fourth rung of the ladder is important, which includes the exact exchange of HF or similar orbital dependent schemes. Another benefit of HF scheme is that it is always free form self interaction error, which is already discussed in Sec. 2. The idea of mixing one part of HF exchange with LDA/GGA/MGGA density functional approximations (DFA) has come from the adiabatic connection formulation (see Sec. 4.4.2 of (engel2013, )) and the exchange-correlation energy in hybrid functional scheme is expressed as:

E_{xc}^{hybrid}=E_{xc}^{DFA}+w~{}\{E_{x}^{HF}-E_{x}^{DFA}\}

(79)

In global hybrids, weight factor $w$ is taken as constant, so, applicable over all space. Some popular global hybrids are B3LYP (B3LYP, ), non-empirical PBE0 (PBE0, ) or TPSSh MGGA (TPSSh, ) hybrids. While global hybrids bring improvement over semi-local DFAs, the inclusion of exact HF type exchange in extended systems is computationally problematic.

To reduce the computational complexity and to achieve better results, the weight factor $w$ \textcolorblackhas been taken as a function of ${\bf r}-{\bf r}^{\prime}$ by separating in short range (SR) and long range (LR) pieces in range separated hybrids (RSH) as:

E_{xc}^{RSH}=E_{xc}^{DFA}+w_{SR}\{E_{x}^{SR-HF}-E_{x}^{SR-DFA}\}+w_{LR}\{E_{x}^{LR-HF}-E_{x}^{LR-DFA}\}

(80)

A popular choice of range separation for the electron-electron repulsion operator is by using general error function $\operatorname{erf}(\vartheta~{}|{\bf r}-{\bf r}^{\prime}|)$ as:

\frac{1}{|{\bf r}-{\bf r}^{\prime}|}=\underbrace{\frac{1-\operatorname{erf}(\vartheta~{}|{\bf r}-{\bf r}^{\prime}|)}{|{\bf r}-{\bf r}^{\prime}|}}_{SR}+\underbrace{\frac{\operatorname{erf}(\vartheta~{}|{\bf r}-{\bf r}^{\prime}|}{|{\bf r}-{\bf r}^{\prime}|}}_{LR}

(81)

Here, the parameter $\vartheta$ controls the range of separation, as, $\vartheta\rightarrow(0/\infty)\equiv(LR/SR)$ . In screened hybrids the LR part is ignored, and, the hybridization with HF like functions is done in short range only. This reduces the computational cost significantly. Heyd, Scuseria and Ernzerhof (HSE) have introduced screened hybrid, which uses the error function based range separation in short range only, and, $w_{SR}^{HSE}=0.25$ (HSE, ). In the modified HSE06 hybrid $\vartheta=0.11~{}a_{0}^{-1}$ , where, $a_{0}$ is the Bohr radius (HSE1, ; HSE06, ).

Short range hybrids are computationally less demanding than the global counterpart for larger systems, but, the short range correction may not provide exact asymptotic nature for molecules. Quantities sensitive to the long-range exchange potential, as well as, to the self-interaction error in density tails may not be addressed well in SR hybrids. Some long range hybrids have been proposed (e.g., LC- $\omega$ PBE (vydrov2007, )) and Henderson, Izmaylov, Scuseria, and Savin (HISS) have proposed a three range separated hybrid, where, a mid range variant is added (HISS, ).

There is yet another genre of hybrids, the local hybrids (LH), where, instead of exchange energies $E_{x}$ , exchange energy densities ( $\tilde{\varepsilon}_{x}$ ) are mixed locally in space:

E_{xc}^{LH}=E_{xc}^{DFA}+\int w({\bf r})~{}\{\tilde{\varepsilon}_{x}^{HF}({\bf r})-\tilde{\varepsilon}_{x}^{DFA}({\bf r})\}d{\bf r}

(82)

The problems with local hybrids, and, the recent developments towards the solution is reviewed in Ref. (maier2019, ), whereas, the different types of hybrids are compared in (henderson2008, ; arbuznikov2007, ). Performance of global and local hybrids in reduction of many electron self-interaction error is compared, and, it has been found that in smaller systems local hybrids perform better. But, with increasing system size, the performance of local hybrid becomes similar to a global hybrids (schmidt2016, ).

5.5 Random Phase Approximation (RPA) within DFT

\textcolor

blackThe random phase approximation scheme comes the fifth rung of the Jacob’s ladder. RPA is a widely used method in different branches of many-body scattering theory. Within DFT scheme, RPA can be formulated using the adiabatic coupling fluctuation-dissipation theorem (ACFDT), though, this is not the only method. RPA provides the long range correlation exactly, so, it is supposed to provide a good description of electronic structures. For molecular systems having complex interaction behaviour, this may come to be handy, but for solids, the necessity of RPA is limited. However, the computational cost for RPA calculation is still so high that for large systems it is not widely used yet. For more mathematical details on ACFDT-RPA method, readers are referred to (ren2012, ).

6 Orbital Based Exchange Correlation

In KS-DFT the ground state energy $E[\rho_{0}]$ is a functional of $\rho_{0}({\bf r})$ , and, the effective potential of $N$ electronic $S$ system $v_{eff}({\bf r})$ (Eq. [47]) is found by variational minimization of $E[\rho]$ (Eq. [30]) with respect to variation of $\rho({\bf r})$ . Since, in $S$ system, spin-orbitals $\{\psi_{i}({\bf x})\}$ are functional of $\rho({\bf r})$ , so, $E$ can also be expressed as a functional of spin-orbitals as $E[\{\psi_{i}\}]$ , and, energy minimization can also be achieved through the variation of $E[\{\psi_{i}\}]$ with respect to $v_{eff}({\bf r})$ . This is the basis of Optimised Potential Method (OPM) (see, Chap-2 of (fiolhais2003, )). In OPM, an integral equation which determines $v_{eff}$ (including effect of $v_{xc}$ ) is solved simultaneously with the KS equation. Once again, let, the spin-orbitals be written as combination of space and spin parts: $\psi_{i}({\bf x})=\phi_{i}({\bf r})\chi_{1}(\sigma)$ .

Because the expression of exchange energy $E_{x}$ for $S$ system is not explicitly known in term of the density, the exchange potential $v_{x}({\bf r})$ can be expressed in term of the orbitals as:

v_{x}([\rho],{\bf r})=\frac{\delta E_{x}[\{\phi_{i}\}]}{\delta\rho({\bf r})}=\sum_{i}^{occ}\iint\left[\frac{\delta E_{x}[\{\phi_{i}\}]}{\delta\phi_{i}({\bf r}^{\prime\prime})}\frac{\delta\phi_{i}({\bf r}^{\prime\prime})}{\delta v_{eff}({\bf r}^{\prime})}+c.c.\right]\frac{\delta v_{eff}({\bf r}^{\prime})}{\delta\rho({\bf r})}d{\bf r}^{\prime\prime}d{\bf r}^{\prime}

(83)

Now, $\frac{1}{\phi_{i}^{\ast}({\bf r})}\frac{\delta E_{x}[\{\phi_{i}\}]}{\delta\phi_{i}({\bf r})}=v_{x,i}({\bf r})$ as in HF scheme (Eq. [19]) and in OPM it is explicitly known. Whereas, $\frac{\delta\phi_{i}({\bf r}^{\prime})}{\delta v_{eff}({\bf r})}=-G_{i}({\bf r}^{\prime},{\bf r})\phi_{i}({\bf r})$ is calculated using first-order perturbation theory, where, $G_{i}({\bf r},{\bf r}^{\prime})=\sum\limits_{j,i\neq j}\frac{\phi_{j}({\bf r})\phi_{j}^{\ast}({\bf r}^{\prime})}{\epsilon_{j}-\epsilon_{i}}$ is the Green’s function, $\{\epsilon_{i}\}$ are orbital energies. Inverse of $\frac{\delta v_{eff}({\bf r}^{\prime})}{\delta\rho({\bf r})}$ is the static response function:

\varUpsilon({\bf r},{\bf r}^{\prime})=-\sum\limits_{i}^{occ}\phi_{i}^{\ast}({\bf r})G_{i}({\bf r},{\bf r}^{\prime})\phi_{i}({\bf r}^{\prime})+c.c.

As mentioned, in OPM the total energy is minimised by varying the effective potential, so that, the minimization condition becomes: $\frac{\delta E[\{\psi_{i}\}]}{\delta v_{eff}({\bf r})}=0$ . As, $v_{ext}({\bf r})$ and $v_{H}({\bf r})$ are exactly known, the minimization condition leads to the integral equation:

\sum_{i}\int[v_{x}^{OPM}({\bf r}^{\prime})-v_{x,i}({\bf r}^{\prime})]\phi_{i}^{\ast}({\bf r}^{\prime})G_{i}({\bf r}^{\prime},{\bf r})\phi_{i}({\bf r})d{\bf r}^{\prime}+c.c.=0

(84)

This integral equation is to be solved to get the minimised energy in OPM.

Kotani has first carried out systematic study on exact exchange potentials on a variety of semiconductors using OPM method and these studies have showed excellent bandgap matching with experiments (kotani1995, ; kotani1996, ).

Another orbital dependent approach comes from the quantal DFT approach (sahni2016, ). In Harbola-Sahni (HS) potential method, the exchange potential is expressed as a line integral over an electric field ${\bf\mathcal{E}}({\bf r})$ (harbola1989, ; sahni1990, ):

W^{HS}({\bf r})=-\int_{\infty}^{{\bf r}}\mathcal{E}({\bf r}^{\prime}).{\bf dl}~{};~{}\text{where,}~{}~{}\mathcal{E}=\int\frac{\rho_{x}({\bf r},{\bf r}^{\prime})}{|{\bf r}-{\bf r}^{\prime}|^{3}}({\bf r}-{\bf r}^{\prime})d{\bf r}^{\prime}

(85)

For spherically symmetric densities, as used in atomic sphere approximation (ASA) in LMTO-ASA, this $\mathcal{E}$ is curl free. For nonspherical charge densities, the solenoidal part of $\mathcal{E}$ is related to the difference in kinetic energies between the HF and the HS approaches (sahni1997, ) and the contribution is numerically insignificant (slamet1994, ). The HS potential can directly be derived from Schrödinger equation (holas1997, ). Using virial theorem, the exchange energy is found as:

E_{x}^{HS}=\int\rho({\bf r})~{}{\bf r}.\nabla W^{HS}({\bf r})d{\bf r}

(86)

7 Performance Comparison of different DFT Methods

DFT methods are meant to reduce the computational costing, yet producing exact physical properties of materials (giustino2014, ). We have seen that, it is a hard task to balance between the atomic-molecular system and solid state systems, as, the nature of the variation of the charge densities are fundamentally different. For molecular systems, a rigorous analysis on the performance of different methods can be found in many texts, e.g, (mardirossian2017, ). Here we focus on the solids, specifically, semiconductors.

To compare the performance of different functionals, we are relying on two parameters, mean percentage error MPE= $\frac{100\%}{N}\sum\limits_{i}^{N}\frac{x_{i}-x_{i}^{expt}}{x_{i}^{expt}}$ , and, mean absolute percentage error MAPE= $\frac{100\%}{N}\sum\limits_{i}^{N}\frac{|x_{i}-x_{i}^{expt}|}{|x_{i}^{expt}|}$ , where, $\{x_{i}\}$ are DFT calculated values. In case of more than one experimental values, the average of those are taken as $x_{i}^{expt}$ . Having large MPE or MAPE means large deviation from experimental values, indicating poor performance. Now, there is another matter of concern, the consistency of results. If the absolute value of MPE differs much from MAPE, then it indicates that the result from this particular functional sometimes give larger than experimental values, and, sometimes smaller than that, thus, reducing the consistency. As an example, we can take the case of PBE GGA; it overestimates the lattice constants, and, largely underestimates the bandgap of semiconductors, although, the results are consistent in the sense that, one almost never find a PBE lattice constant lower than the experimental value and a PBE bandgap larger than the experimental value.

Another important thing to mention here is that the implementation of same functional in different DFT packages may produce different results, and in some cases, the implementation is tricky, and, sometimes become package dependent as well. A recent article is on the reproducibility of DFT results has aimed to target this issue, and, worth reading (lejaeghere2016, ), although, discussion on every single package can not be included in a paper anyhow. So, we try to mention the method and/or package in the tables wherever possible.

7.1 Structural Parameters

Lattice Constant

Optimized structure is the basic criteria for further study on any material. Geometric optimization is done using energy minimization technique, where variation of energy is found by varying the lattice parameters of any solid (variable cell relaxation), as well as, by changing relative positions of atoms (ionic relaxation). Thus, the energy minimised structure seems to be most stable structure. Although LDA is a very basic approximation, it can predict the lattice constants quite well. This is because the HEG or its slowly varying approximation is suitable for solids, though for molecules it is not realistic. PBE GGA incorporates the slowly varying density idea to LDA, so, performs better than LDA (see Table 7.1). One can see from this table, MPE and MAPE are same for PBE, so, PBE consistently overestimates the lattice constants for semiconductor. Now, let us take the case of Ge. The deviation is quite high, i.e., PBE over-binds the atoms so much. This is because of the deviation of PBE from GEA result, as we have discussed. PBEsol is better suited for solids, and, the MPE and MAPE proves this simple correction is one of the best performer in the test on lattice constants, even in this modern generation of DFT functionals. First non-empirical MGGA, TPSS is built upon PBE, so, overestimation is also evident. However, the SCAN MGGA functional can reduce the error by appropriate constraint management, but, the consistency as seen in PBE or TPSS, is not that much robust in SCAN. TM reduces the error significantly, and, the consistency is manageable (jana2018, ; mo2017, ). In Fig. 6 the MPE for different semiconductors are plotted taking the values from (jana2018, ).

	$\displaystyle E^{HF}_{ee}=\langle\Phi^{HF}\|\hat{V}_{ee}\|\Phi^{HF}\rangle$
$\displaystyle=\iint\frac{d{\bf r}_{1}d{\bf r}_{2}}{2\|{\bf r}_{1}-{\bf r}_{2}\|}$	$\displaystyle\left[\phi^{\ast}_{1}({\bf r}_{1})\phi^{\ast}_{2}({\bf r}_{2})\{\phi_{1}({\bf r}_{1})\phi_{2}({\bf r}_{2})-\phi_{2}({\bf r}_{1})\phi_{1}({\bf r}_{2})\}\right.$
	$\displaystyle\left.-\phi^{\ast}_{2}({\bf r}_{1})\phi^{\ast}_{1}({\bf r}_{2})\{\phi_{1}({\bf r}_{1})\phi_{2}({\bf r}_{2})-\phi_{2}({\bf r}_{1})\phi_{1}({\bf r}_{2})\}\right]$
$\displaystyle=\iint\frac{d{\bf r}_{1}d{\bf r}_{2}}{2\|{\bf r}_{1}-{\bf r}_{2}\|}$	$\displaystyle\{\|\phi_{1}({\bf r}_{1})\|^{2}\|\phi_{2}({\bf r}_{2})\|^{2}+\|\phi_{2}({\bf r}_{1})\|^{2}\|\phi_{1}({\bf r}_{2})\|^{2}\}$
$\displaystyle-\iint\frac{d{\bf r}_{1}d{\bf r}_{2}}{2\|{\bf r}_{1}-{\bf r}_{2}\|}$	$\displaystyle\{\phi^{\ast}_{1}({\bf r}_{1})\phi_{1}({\bf r}_{2})\phi^{\ast}_{2}({\bf r}_{2})\phi_{2}({\bf r}_{1})-\phi^{\ast}_{1}({\bf r}_{2})\phi_{1}({\bf r}_{1})\phi^{\ast}_{2}({\bf r}_{1})\phi_{2}({\bf r}_{2})\}$	(4)

	$\displaystyle E_{ee}^{HF}=\langle\Psi^{HF}\|\hat{V}_{ee}\|\Psi^{HF}\rangle=\frac{1}{2}$	$\displaystyle\sum_{i,j=1}^{N}\left(J_{ij}-K_{ij}\right)=E_{H}+E_{x}$		(8)
	$\displaystyle J_{ij}=\iint\frac{\psi_{i}^{\ast}({\bf x})\psi_{i}({\bf x})\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{j}({\bf x}^{\prime})}{\|{\bf r}-{\bf r}^{\prime}\|}d{\bf x}d{\bf x}^{\prime}$	$\displaystyle;K_{ij}=\iint\frac{\psi_{i}^{\ast}({\bf x})\psi_{i}({\bf x}^{\prime})\psi_{j}^{\ast}({\bf x}^{\prime})\psi_{j}({\bf x})}{\|{\bf r}-{\bf r}^{\prime}\|}d{\bf x}d{\bf x}^{\prime}$

	$\displaystyle\langle\Psi^{HF}\|\hat{T}+\hat{V}_{ext}\|\Psi^{HF}\rangle$	$\displaystyle+2E_{ee}^{HF}=\sum_{i}^{N}\epsilon_{i}$
	$\displaystyle\Rightarrow E[\Psi^{HF}]=\langle\Psi^{HF}\|\hat{T}$	$\displaystyle+\hat{V}_{ext}\|\Psi^{HF}\rangle+E_{ee}^{HF}=\sum_{i}^{N}\epsilon_{i}-E_{ee}^{HF}$		(18)

		$\displaystyle\text{KS: }~{}~{}~{}~{}~{}~{}F^{S}[\rho_{0}]=\inf_{\Psi^{t}\rightarrow\rho_{0}({\bf r})}\langle\Psi^{t}\|\hat{T}\|\Psi^{t}\rangle$		(50)
		$\displaystyle\text{HF KS: }~{}~{}~{}~{}~{}~{}F^{S}[\rho_{0}]=\inf_{\Psi^{t}\rightarrow\rho_{0}({\bf r})}\left(\langle\Psi^{t}\|\hat{T}\|\Psi^{t}\rangle+E_{H}[\{\psi_{i}\}]+E_{x}[\{\psi_{i}\}]\right)$
		$\displaystyle\text{sx-LDA: }~{}~{}~{}~{}F^{S}[\rho_{0}]=\inf_{\Psi^{t}\rightarrow\rho_{0}({\bf r})}\left(\langle\Psi^{t}\|\hat{T}\|\Psi^{t}\rangle+E_{H}[\{\psi_{i}\}]+E^{sx}_{x}[\{\psi_{i}\}]\right)$

Methods $\rightarrow$ Solids $\downarrow$	LSDA	PBE	PBEsol	TPSS	SCAN	TM	Expt.
\svhline C (A4)	3.536	3.573	3.557	3.572	3.555	3.554	3.567
Si (A4)	5.400	5.467	5.433	5.450	5.425	5.411	5.430
Ge (A4)	5.648	5.785	5.704	5.754	5.687	5.672	5.652
SiC (B3)	4.332	4.379	4.359	4.365	4.352	4.344	4.358
BN (B3)	3.583	3.625	3.607	3.624	3.605	3.608	3.607
BP (B3)	4.490	4.546	4.521	4.545	4.521	4.510	4.538
BAs (B3)	4.742	4.817	4.778	4.810	4.779	4.763	4.777
AlP (B3)	5.433	5.504	5.470	5.489	5.478	5.450	5.460
AlAs (B3)	5.637	5.732	5.681	5.707	5.670	5.656	5.658
AlSb (B3)	6.120	6.232	6.168	6.208	6.173	6.143	6.136
$\beta-$ GaN (B3)	4.503	4.588	4.547	4.581	4.524	4.549	4.531
GaP (B3)	5.425	5.533	5.474	5.523	5.457	5.464	5.448
GaAs (B3)	5.627	5.763	5.684	5.737	5.664	5.664	5.648
GaSb (B3)	6.067	6.226	6.130	6.190	6.117	6.102	6.096
InP (B3)	5.878	6.001	5.932	5.989	5.938	5.923	5.866
InAs (B3)	6.061	6.211	6.122	6.182	6.122	6.104	6.054
InSb (B3)	6.472	6.651	6.543	6.611	6.545	6.521	6.479
ZnS (B3)	5.403	5.440	5.355	5.401	5.370	5.364	5.409
ZnSe (B3)	5.570	5.734	5.634	5.681	5.652	5.633	5.668
ZnTe (B3)	5.995	6.178	6.064	6.115	6.077	6.056	6.089
CdS (B3)	5.758	5.926	5.824	5.933	5.856	5.857	5.818
CdSe (B3)	6.009	6.195	6.080	6.192	6.100	6.102	6.052
CdTe (B3)	6.405	6.610	6.291	6.604	6.521	6.497	6.480
MgO (B1)	4.145	4.242	4.206	4.224	4.184	4.202	4.207
MgS (B3)	5.580	5.684	5.642	5.681	5.634	5.629	5.202
MgSe (B1)	5.382	5.501	5.445	5.491	5.454	5.435	5.400
MgTe (B3)	6.365	6.506	6.439	6.500	6.452	6.422	6.420
CaS (B1)	5.570	5.710	5.632	5.698	5.683	5.657	5.689
CaSe (B1)	5.798	5.955	5.869	5.947	5.921	5.894	5.916
CaTe (B1)	6.215	6.389	6.291	6.386	6.375	6.317	6.348
SrS (B1)	5.910	6.056	5.973	6.047	6.031	6.007	5.990
SrSe (B1)	6.129	6.297	6.203	6.286	6.264	6.234	6.234
SrTe (B1)	6.531	6.714	6.609	6.708	6.693	6.641	6.640
BaS (B1)	6.289	6.433	6.362	6.448	6.441	6.390	6.389
BaSe (B1)	6.510	6.681	6.577	6.670	6.659	6.622	6.595
BaTe (B1)	6.890	7.080	6.964	7.075	7.071	7.012	7.007
MPE (% )	-0.678	1.502	0.186	1.238	0.556	0.271	–
MAPE (% )	1.099	1.502	0.789	1.247	0.723	0.588	–

Semiconductor Physics: A Density Functional Journey

Abstract

1 Introduction

2 The Hartree - Fock Method

2.1 Hartree - Fock Method for $N$ electron system

2.2 Hartree - Fock Equation

2.3 Hartree - Fock - Slater Equation

3 Hohenberg-Kohn (HK) Theorem

4 The Kohn-Sham (KS) DFT

4.1 Generalised KS Schemes

4.2 KS DFT for Fractional $N$ ; Derivative Discontinuity and Bandgap:

5 Jacob’s Ladder of Exchange-Correlation Functional

5.1 Local Density Approximation (LDA)

5.2 Generalised Gradiant Approximation (GGA)

5.3 Meta-GGA

5.4 Hybrid Functional Method

5.5 Random Phase Approximation (RPA) within DFT

6 Orbital Based Exchange Correlation

7 Performance Comparison of different DFT Methods

7.1 Structural Parameters

Lattice Constant

Bulk Modulus

Transition Pressure for Structural Phase Transition

References

$\hat{H}\rightarrow$ Hamiltonion Operator	$\hat{T}$ & $\hat{V}\rightarrow$ Kin. & Pot. energy Operator
$\Psi({\bf X})\rightarrow$ Spin dependent wavefn.	$\Phi({\bf R})\rightarrow$ Spin independent wavefn.
$\Psi^{S}({\bf X})\rightarrow$ Slater ,,	$\Phi^{S}({\bf R})\rightarrow$ Slater ,,
$\psi_{i}({\bf x}_{n})=\phi_{i}({\bf r}_{n})\chi_{i}(\sigma_{n})\rightarrow$ $i$ -th Spin-orbital state at ${\bf r}_{n}$ and of spin $\sigma_{n}$	$\phi_{i}({\bf r}_{n}),\chi_{i}(\sigma_{n})\rightarrow$ $i$ -th Space-only orbital at ${\bf r}_{n}$ and $i$ -th spin state.
$\rho({\bf r}),\rho_{0}({\bf r})\rightarrow$ Density and Ground State Density	$\rho_{1}({\bf r},{\bf r}^{\prime}),\rho_{2}({\bf r}_{1}{\bf r}_{2},{\bf r}_{1}^{\prime}{\bf r}_{2}^{\prime})\rightarrow$ 1st and 2nd order Spinless Reduced Density Matrix.
$\rho_{x}({\bf r},{\bf r}^{\prime})\rightarrow$ Fermi Exchange Hole	$\rho_{xc}({\bf r},{\bf r}^{\prime})\rightarrow$ Fermi-Coulomb xc Hole
$g_{xc}\rightarrow$ Pair-Correlation Function	$\mu_{\pm}\rightarrow$ Chemical Potential
$I^{(N)}\rightarrow$ Ionization pot. of $N$ $e$ system.	$A^{(N)}\rightarrow$ Electron Affinity.
$E\rightarrow$ Total Energy; $E_{g}\rightarrow$ Bandgap	$E_{ee}\rightarrow$ e-e Interaction Energy
$E_{H}\rightarrow$ Hartree Energy	$E_{xc}\rightarrow$ Exchange-Correlation Energy
$\epsilon_{i}\rightarrow$ $i$ -th orbital energy	$\varepsilon_{xc}\rightarrow$ xc Energy per particle.
$\tau({\bf r})\rightarrow$ Kinetic Energy Density	$\tilde{f}\rightarrow$ Exchange Enhancement Factor.
$v_{eff},v_{ext}\rightarrow$ Effective and External Pot.	$v_{H},v_{xc}\rightarrow$ Hartree and xc Pot.

Semiconductor Physics: A Density Functional Journey

Abstract

1 Introduction

2 The Hartree - Fock Method

2.1 Hartree - Fock Method for NN electron system

2.2 Hartree - Fock Equation

2.3 Hartree - Fock - Slater Equation

3 Hohenberg-Kohn (HK) Theorem

4 The Kohn-Sham (KS) DFT

4.1 Generalised KS Schemes

4.2 KS DFT for Fractional NN; Derivative Discontinuity and Bandgap:

5 Jacob’s Ladder of Exchange-Correlation Functional

5.1 Local Density Approximation (LDA)

5.2 Generalised Gradiant Approximation (GGA)

5.3 Meta-GGA

5.4 Hybrid Functional Method

5.5 Random Phase Approximation (RPA) within DFT

6 Orbital Based Exchange Correlation

7 Performance Comparison of different DFT Methods

7.1 Structural Parameters

Lattice Constant

Bulk Modulus

Transition Pressure for Structural Phase Transition

References

2.1 Hartree - Fock Method for $N$ electron system

4.2 KS DFT for Fractional $N$ ; Derivative Discontinuity and Bandgap: