\chapter{Theoretical background: density functional theory} \label{ch:chapter3}
\lhead{\emph{Theoretical background: density functional theory}}

In this Chapter we outline the theoretical methods used in our analysis of the structural, elastic and electronic properties of the III-P compounds AlP, GaP and InP in Chapter~\ref{ch:chapter6}. We begin in Sec.~\ref{sec:ch_3_kohn_sham} by describing the Kohn-Sham formalism of density functional theory (DFT). In Sec.~\ref{sec:ch_3_impl}, we provide an outline of key considerations and approximations employed in the implementation of Kohn-Sham DFT calculations. In Sec.~\ref{sec:ch_3_impl_LDA_TB-mBJ}, we begin by discussing schemes for the determination of the exchange-correlation (XC) energy in Kohn-Sham DFT, such as the local density approximation (LDA), generalised gradient approximations (GGAs) and the Tran-Blaha modified Becke-Johnson (TB-mBJ) exchange potential. In Sec.~\ref{sec:ch_3_impl_HSE}, we outline another approach for the determination of the XC energy, the hybrid Heyd-Scuseria-Ernzerhof functional for solids (HSEsol). In Sec.~\ref{sec:ch_3_impl_PAW}, we overview the use of pseudopotentials for the description of the Kohn-Sham effective potential, and the projector augmented-wave (PAW) method for the description of all-electron wave functions. In Sec.~\ref{sec:ch_3_impl_plane_waves}, we describe the plane wave formalism for Kohn-Sham DFT calculations. Finally, in Sec.~\ref{sec:ch_3_conclusions} we summarise and conclude.

We note that our analysis of III-P band parameters by means of DFT calculations in Chapter~\ref{ch:chapter6}, using the methods outlined in this Chapter, are implemented using the Vienna Ab-initio Software Package (VASP) \cite{vasp_website,Kresse_PRB_1996}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Section 1: Kohn-Sham DFT %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Kohn-Sham density functional theory}
\label{sec:ch_3_kohn_sham}

% Kohn-Sham

In this Section we present an outline of the Kohn-Sham formulation of DFT, comprehensive detail of which is provided by Martin \cite{Martin_book_2004}. The Hohenberg-Kohn theorems outlined that the electron density $n ( \textbf{r} )$ can be considered as a ``basic variable'', such that any system of many interacting particles can be viewed as a functional of the ground state density $n_{0}(\textbf{r})$ \cite{Hohenberg_PR_1964}. As such, this density determines all the information in the many-body wave function of the ground state. Kohn and Sham proposed an independent-particle approach for the determination of the density \cite{Kohn_PR_1965}. This approach, building upon the Hohenberg-Kohn theorems, approximates the density $n(\textbf{r})$ via a set of independent wave functions $\psi_{n} ( \textbf{r} ) = \langle \textbf{r} \vert \psi_{n} \rangle$ as

\begin{equation}
    n (\textbf{r}) = \sum_{n \in \text{occ}} \vert \psi_{n} (\textbf{r}) \vert^{2} \, , \label{eq:ch_3_KS_charge_density}
\end{equation}

\noindent
where $n$ runs over all occupied states. In the Kohn-Sham formalism, all many-body effects are assumed to be encapsulated by an XC functional, however no exact functionals are known for any system of more than one electron \cite{Martin_book_2004}. We now outline the evaluation of the total energy of the a many-body system using the Kohn-Sham formalism.

% Total energy

In calculating the electronic ground state, the total energy of the system must be evaluated. The Hamiltonian operator for a system consisting of an arbitrary number of interacting electrons and nuclei can be written as \cite{Martin_book_2004}

\begin{equation}
    \widehat{H} = \widehat{T} + \widehat{V}_{\text{ext}} + \widehat{V}_{\text{int}} + \widehat{T}_{\text{nuc}} + E_{\text{II}} \, , \label{eq:ch_3_KS_hamiltonian_martin}
\end{equation}

\noindent
where $\widehat{T}$ is the kinetic operator for the electrons, $\widehat{V}_{\text{ext}}$ is the attractive electrostratic potential acting on the electrons due to the nuclei, $\widehat{V}_{\text{int}}$ describes the repulsive electron-electron electrostatic interactions, $\widehat{T}_{\text{nuc}}$ is the kinetic operator for the nuclei, and $E_{\text{II}}$ describes the repulsive nucleus-nucleus interactions. The total energy of the system is given by the expectation value of Eq.~\eqref{eq:ch_3_KS_hamiltonian_martin}

\begin{equation}
    E = \bra{\psi} \widehat{H} \ket{\psi} \equiv \langle \widehat{H} \rangle \, , \label{eq:ch_3_total_energy_expectation_value}
\end{equation}

\noindent
where $\psi ( \textbf{r} ) = \langle \textbf{r} \vert \psi \rangle$ is the normalised wave function characterising the state of the many-electron system.

When providing an explicit form of Eq.~\eqref{eq:ch_3_KS_hamiltonian_martin}, approximations can be made in order to simplify computation \cite{Martin_book_2004,Caro_thesis_2013}. The first of these is the Born-Oppenheimer approximation \cite{Born_AdP_1927}. This approximation considers the movement of electrons and nuclei separately, and assumes that the nuclear positions are fixed. This approximation is made as the more massive nuclei move slowly compared with the motion of electrons. This allows for the separation of the wave function into parts containing the nuclear and electronic degrees of freedom independently as

\begin{equation}
    \ket{\psi} = \ket{\psi}_{e} \times \ket{\psi}_{I} \, , \label{eq:ch_3_wave_function_born_opp_approx}
\end{equation}

\noindent
where $\ket{\psi}_{e}$ $\left( \ket{\psi}_{I} \right)$ is the electron (nuclei) wavefunction. In the Born-Oppenheimer approximation, we focus on $\ket{\psi}_{e}$, which for simplicity we denote as $\ket{\psi}_{e} \equiv \ket{\psi}$. As such, in Eq.~\eqref{eq:ch_3_KS_hamiltonian_martin} we neglect $\widehat{T}_{\text{nuc}}$, the nuclei kinetic operator. This imposes that the nucleus-nucleus interactions $E_{II}$ become constant \cite{Caro_thesis_2013} -- i.e.~ simply a reference energy -- with the nuclear positions becoming parameters. The total energy can then be written as

\begin{equation}
    E = \langle \widehat{T} \rangle + \int V_{\text{ext}} (\textbf{r}) \, n (\textbf{r}) \, \text{d}\textbf{r} + \langle \widehat{V}_{\text{int}} \rangle \, . \label{eq:ch_3_total_energy_approx}
\end{equation}

The calculation of the expectation value of $\widehat{V}_{\text{int}}$, $\langle \widehat{V}_{\text{int}} \rangle$, is one of the main difficulties in DFT calculations \cite{Caro_thesis_2013}. As an approximation, $\langle \widehat{V}_{\text{int}} \rangle$ can be expressed as the sum of the classical electrostatic energy and the XC energy $E_{\scalebox{0.6}{\text{XC}}}$

\begin{equation}
    \langle \widehat{V}_{\text{int}} \rangle = E_{\scalebox{0.6}{\text{XC}}} + \frac{1}{8\pi \epsilon_{0}} \int \int \frac{ n({\textbf{r}}) n(\textbf{r}') }{ \vert \textbf{r} - \textbf{r}' \vert } \, \text{d}\textbf{r} \, \text{d}\textbf{r}' \, , \label{eq:ch_3_expectation_V_internal}
\end{equation}

\noindent
where $E_{\scalebox{0.6}{\text{XC}}}$ is the difference between the energy of the quantum-mechanical electron-electron interaction and that corresponding to the classical electrostatic electron-electron interaction, or Hartree energy $E_{\scalebox{0.6}{\text{Hartree}}}$ \cite{Caro_thesis_2013}. Initial attempts to solve Eq.~\eqref{eq:ch_3_total_energy_approx} by Thomas \cite{Thomas_MPCPS_1927} and Fermi \cite{Fermi_RANL_1927}, leading to the Thomas-Fermi model, produced large errors compared to calculations that neglect $E_{\scalebox{0.6}{\text{XC}}}$ \cite{Caro_thesis_2013}. To overcome the limitations of the Thomas-Fermi model, Kohn and Sham presented their formalism, as described in Eq.~\eqref{eq:ch_3_KS_charge_density}, where the density is constructed via a set of independent single-electron wave functions (the so-called Kohn-Sham orbitals). Using Eq.~\eqref{eq:ch_3_KS_charge_density}, the electron kinetic energy can be expressed as

\begin{equation}
    T_{\scalebox{0.6}{\text{KS}}} = - \frac{\hbar^{2}}{2m_{0}} \sum_{n \in \text{occ}} \bra{\psi_{n}} \nabla_{n}^{2} \ket{\psi_{n}} \, , \label{eq:ch_3_KS_kinetic_energy}   
\end{equation}

\noindent
so that the total energy can be written as

\begin{eqnarray}
    E_{\scalebox{0.6}{\text{KS}}} = &-& \frac{\hbar^{2}}{2m_{0}} \sum_{n \in \text{occ}} \bra{\psi_{n}} \nabla_{n}^{2} \ket{\psi_{n}} + \int V_{\text{ext}}(\textbf{r}) \, n(\textbf{r}) \, \text{d}\textbf{r} + \frac{1}{8\pi \epsilon_{0}} \int \int \frac{ n({\textbf{r}}) n(\textbf{r}') }{ \vert \textbf{r} - \textbf{r}' \vert } \, \text{d}\textbf{r} \, \text{d}\textbf{r}' \nonumber \\ &+& E_{\scalebox{0.6}{\text{XC}}} [ n(\textbf{r}) ] + E_{\text{II}} \, . \label{eq:ch_3_KS_total_energy}
\end{eqnarray}

Having now expressed the total energy in the Kohn-Sham formalism, what remains is to determine the XC energy. Following from above, the XC energy can be written as

\begin{equation}
    E_{\scalebox{0.6}{\text{XC}}} = \langle \widehat{T} \rangle - T_{\scalebox{0.6}{KS}} + \langle \widehat{V}_{\text{int}} \rangle - E_{\scalebox{0.6}{\text{Hartree}}} \, , \label{eq:ch_3_KS_XC_functional}
\end{equation}

\noindent
which expresses $E_{\scalebox{0.6}{\text{XC}}}$ as the difference between the energy of the exact many-body interacting system and its energy calculated in the independent-electron approximation. Under application of the variational principle and treating the energy as a functional of the electron density, $E [ n ( \textbf{r} ) ]$, the Kohn-Sham eigenvalue equation, Hamiltonian $\widehat{H}_{\scalebox{0.6}{\text{KS}}}$, and effective potential $\widehat{V}_{\scalebox{0.6}{\text{KS}}}$ are given as \cite{Martin_book_2004}

\begin{equation}
    \widehat{H}_{\scalebox{0.6}{\text{KS}}} \psi_{i} (\textbf{r}) = \epsilon_{i} \psi_{i} \, , \label{eq:ch_3_KS_eigenevalue_eqn}
\end{equation}

\begin{equation}
    \widehat{H}_{\scalebox{0.6}{\text{KS}}} = - \frac{\hbar^{2}}{2m_{0}} \nabla^{2} + \widehat{V}_{\scalebox{0.6}{\text{KS}}} (\textbf{r}) \, , \label{eq:ch_3_KS_hamiltonian}
\end{equation}

\begin{equation}
    \widehat{V}_{\scalebox{0.6}{\text{KS}}} = \widehat{V}_{\text{ext}} (\textbf{r}) + \frac{\delta E_{\scalebox{0.6}{\text{Hartree}}}}{\delta n (\textbf{r})} + \frac{\delta E_{\scalebox{0.6}{\text{XC}}}}{\delta n (\textbf{r})} \, , \label{eq:ch_3_KS_effective_potential}
\end{equation}

\noindent
where $\delta$ denotes the functional derivative. From Eq.~\eqref{eq:ch_3_KS_XC_functional}, we see that the accuracy of the electronic ground state obtained from the implementation of the Kohn-Sham formulation of DFT relies heavily on the choice of XC functional $E_{\scalebox{0.6}{\text{XC}}}  [ n ( \textbf{r} ) ]$, the choice of which is non-trivial. Information regarding the choice of XC functionals employed in this work are provided in Sec.~\ref{sec:ch_3_impl}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Section 2: Pseudopotentials and PAW method %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Considerations in the implementation of Kohn-Sham density functional theory}
\label{sec:ch_3_impl}

In this Section we provide a basic outline of the considerations and techniques employed in our implementation of the Kohn-Sham formalism of DFT, outlined in Sec.~\ref{sec:ch_3_kohn_sham}. We begin in Sec.~\ref{sec:ch_3_impl_LDA_TB-mBJ} by describing the LDA. In Sec.~\ref{sec:ch_3_impl_LDA_TB-mBJ} we also outline the TB-mBJ exchange potential, which with the LDA correlation can be combined to perform DFT calculations in the meta-GGA. Next, in Sec.~\ref{sec:ch_3_impl_HSE} we detail hybrid functionals, specifically the HSEsol hybrid functional, as a generalisation of the Kohn-Sham approach. Following this, in Sec.~\ref{sec:ch_3_impl_PAW} we describe the PAW method, which is a robust approach allowing for the computation of the many-electron ground state via DFT. Finally, in Sec.~\ref{sec:ch_3_impl_plane_waves} we outline the plane wave formalism for implementing DFT calculations.

\subsection{The local density approximation and Tran-Blaha modified Becke-Johnson exchange potential}
\label{sec:ch_3_impl_LDA_TB-mBJ}

% LDA

As discussed in Sec.~\ref{sec:ch_3_kohn_sham}, the choice of XC functional is key in determining the accuracy of Kohn-Sham DFT calculations. One approach for determining the XC energy is the LDA, which builds on analysis of a uniform electron gas (UEG), where the system is completely specified by a constant density \cite{Perdew_PRB_1992,Martin_book_2004,Caro_thesis_2013}

\begin{equation}
    n = \frac{3}{4 \pi r_{s}^{3}} \, , \label{eq:ch_3_ueg_density}
\end{equation}

\noindent
where $r_{s}$ is the Wigner-Seitz radius which defines a sphere containing one electron on average. Thus $r_{s}$ measures the average distance between electrons in a UEG. In the LDA the XC energy is treated as an integral of the product of the electron density and the XC energy density over all space. In the LDA the XC energy density at each point is assumed to be equal to that in a UEG with density equal to the local density, $n(\textbf{r})$. The LDA exchange energy $E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{LDA}}}$ is then given as a functional of the local charge density $n(\textbf{r})$ as \cite{Martin_book_2004,Caro_thesis_2013}

% The LDA is a reduced form of the local spin density approximation, where the system is unpolarised and the spin dependent charge densities are simplified setting $n^{\uparrow}(\textbf{r}) = n^{\downarrow}(\textbf{r}) = n(\textbf{r}) / 2$.

\begin{equation}
    E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{LDA}}} \left[ n(\textbf{r}) \right] = \int n(\textbf{r}) \epsilon_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{UEG}}} \left[ n(\textbf{r}) \right] \text{d}\textbf{r} \, , \label{eq:ch_3_E_X_LDA}
\end{equation}

\noindent
where $\epsilon_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{UEG}}} = E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{UEG}}} / N_{e}$ is the exchange energy of the UEG per electron, with $N_{e}$ being the total number of electrons, and $\epsilon_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{UEG}}} \propto r_{s}^{-1}$. The remaining aspect of the LDA implementation is then to parametrise $\epsilon_{\scalebox{0.6}{\text{XC}}}^{\scalebox{0.6}{\text{UEG}}}$, which can be separated into exchange and correlation contributions, $\epsilon_{\scalebox{0.6}{\text{XC}}}^{\scalebox{0.6}{\text{UEG}}} = \epsilon_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{UEG}}} + \epsilon_{\scalebox{0.6}{\text{C}}}^{\scalebox{0.6}{\text{UEG}}}$. The exchange contribution $\epsilon_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{UEG}}}$ takes an analytical form \cite{Martin_book_2004,Caro_thesis_2013}, while the correlation contribution $\epsilon_{\scalebox{0.6}{\text{C}}}^{\scalebox{0.6}{\text{UEG}}}$ must be numerically parametrised \cite{Martin_book_2004,Caro_thesis_2013}. There exist several approaches to parametrise $\epsilon_{\scalebox{0.6}{\text{C}}}^{\scalebox{0.6}{\text{UEG}}}$, one key example being that of Perdew and Zunger, where they fit to quantum Monte Carlo data for the energy of the UEG \cite{Perdew_PRB_1981}.

% GGA, meta-GGA and TB-mBJ

A key limitation of LDA DFT calculations is that they significantly underestimate the band gap of semiconductors and insulators \cite{Tran_PRL_2009}. While the approach of Perdew and Zunger \cite{Perdew_PRB_1981} gives the exact XC energy for an idealised UEG with constant charge density, in real systems the charge density is position dependent \cite{Caro_thesis_2013}, where $\epsilon_{\scalebox{0.6}{\text{XC}}} (\textbf{r})$ is dependent not only on the local density $n ( \textbf{r} )$, but also semi-local contributions encapsulated by $\nabla n ( \text{r} )$ and higher-order spatial derivatives. As such, we now examine methods beyond the LDA for the determination of the XC energy, which allow to address the band gap underestimation problem. In the GGA it is assumed that $\epsilon_{\scalebox{0.6}{\text{XC}}} (\textbf{r})$ is a functional of $n(\textbf{r})$ and $\nabla n (\textbf{r})$, which gives \cite{Perdew_PRL_1997,Martin_book_2004,Caro_thesis_2013}

\begin{equation}
    E_{\scalebox{0.6}{\text{XC}}}^{\scalebox{0.6}{\text{GGA}}} \left[ n (\textbf{r}) \right] = \int n(\textbf{r}) \, \epsilon_{\scalebox{0.6}{\text{XC}}}^{\scalebox{0.6}{\text{GGA}}} \left[n(\textbf{r}) , \nabla n (\textbf{r}) \right] \, \text{d}\textbf{r} \, . \label{eq:ch_3_E_XC_GGA}
\end{equation}

There exist a number of GGAs, such as the parametrisation of Perdew, Burke and Ernzerhof, the PBE functional \cite{Perdew_PRL_1997}. Another approach is the meta-GGA, where an appropriate exchange potential is employed which takes into account the local kinetic energy density $\tau ( \textbf{r} ) \propto \nabla^{2} n ( \textbf{r} )$. In the meta-GGA the XC functional is therefore given by

\begin{equation}
    E_{\scalebox{0.6}{\text{XC}}}^{\scalebox{0.6}{\text{mGGA}}} \left[ n (\textbf{r}) \right] = \int n(\textbf{r}) \, \epsilon_{\scalebox{0.6}{\text{XC}}}^{\scalebox{0.6}{\text{mGGA}}} \left[n(\textbf{r}) , \nabla n (\textbf{r}) , \nabla^{2} n ( \textbf{r} ) \right] \, \text{d}\textbf{r} \, . \label{eq:ch_3_E_XC_meta-GGA}
\end{equation}

The TB-mBJ exchange potential proposed by Tran and Blaha \cite{Tran_PRL_2009}, building on the work of Becke and Johnson \cite{Becke_JCP_2006}, is an example of a meta-GGA approach to DFT. The TB-mBJ exchange potential yields calculated band gaps that achieve accuracy comparable to more computationally expensive approaches \cite{Tran_PRL_2009}, making it an increasingly popular approach to analyse the electronic structure of semiconductor materials. The TB-mBJ exchange potential is \cite{Tran_PRL_2009}

\begin{equation}
    V_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{mBJ}}} (\textbf{r}) = c \, V_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{BR}}} (\textbf{r}) + (3c-2) \frac{1}{\pi} \sqrt{\frac{5}{12}} \sqrt{\frac{2 \tau (\textbf{r})}{n(\textbf{r})}} \, , \label{eq:ch_3_TB_mBJ_exchange_potential}
\end{equation}

\noindent
where $V^{\scalebox{0.6}{\text{BR}}} (\textbf{r})$ is the Becke-Roussel (BR) potential \cite{Becke_PRA_1989}

\begin{equation}
    V^{\scalebox{0.6}{\text{BR}}} (\textbf{r}) = -\frac{1}{b (\textbf{r})} \left( 1 - e^{-x(\textbf{r})} - \frac{1}{2} x(\textbf{r}) e^{-x(\textbf{r})} \right) \, , \label{eq:ch_3_BR_potential}
\end{equation}

\noindent
where $b(\textbf{r}) = \left[ x(\textbf{r}) e^{-x(\textbf{r})} / \left( 8 \pi n(\textbf{r}) \right)  \right]^{1/3}$. The function $x(\textbf{r})$ is determined via the solutions of a position dependent non-linear equation which determines the exchange hole, exchange energy and corresponding exchange energy density \cite{Becke_PRA_1989,Tran_PRL_2009}. $x(\textbf{r})$ is determined via $\tau(\textbf{r})$, $n(\textbf{r})$, $\nabla n(\textbf{r})$, and $\nabla^{2} n(\textbf{r})$. The inclusion of this BR potential is weighted by the BR ``mixing'' parameter $c$, which Tran and Blaha chose to depend linearly on the square root of the average of $\vert \nabla n (\textbf{r}) \vert / n (\textbf{r})$ as \cite{Tran_PRL_2009}

\begin{equation}
    c = \alpha + \beta \left( \frac{1}{\Omega} \int_{\Omega} \frac{\vert \nabla n (\textbf{r}') \vert}{ n (\textbf{r}') } \text{d}\textbf{r}' \right)^{1/2} \, , \label{eq:ch_2_c_parameter}
\end{equation}

\noindent
where $\alpha$ and $\beta$ are adjustable parameters, and $\Omega$ is the unit cell volume. We note that $c$ can be treated as an empirical parameter through which to fit to an experimental band gap, an approach we will employ and analyse in Chapter~\ref{ch:chapter6}. Tran and Blaha included the BR potential in their modifications to the Becke-Johnson potential, in order to model the Coulomb potential generated by the exchange hole, replacing the Slater potential proposed for this purpose by Becke and Johnson \cite{Becke_JCP_2006}. The TB-mBJ exchange potential recovers the LDA for a constant electron density and effectively mimics the behaviour of orbital-dependent potentials \cite{Tran_PRL_2009}. Therefore, the use of the TB-mBJ exchange potential, which is employed in conjunction with the LDA correlation, is a computationally inexpensive approach for determining the XC potential in Kohn-Sham DFT. We note that the TB-mBJ is a ``potential only'' functional, meaning it has no associated XC energy, and is hence not self-consistent with respect to the total energy. The use of the TB-mBJ exchange potential improves the description of band gaps vs.~the LDA alone \cite{Tran_PRL_2009}. However, a key issue with the TB-mBJ XC potential is its tendency to overestimate band edge effective masses \cite{Kim_PRB_2010,Laurien_PRB_2022}. As a result, to describe the electronic structure with even greater accuracy, below we will examine hybrid functional approaches for the determination of the XC energy. Finally, we note that since TB-mBJ is a ``potential only'' functional, it is not possible to employ TB-mBJ to compute Hellmann-Feynmann ionic forces, thus meaning that it is not possible to perform structural relaxation using this method. In Chapter~\ref{ch:chapter6}, we therefore combine LDA calculations of structural properties with TB-mBJ calculations to analyse the electronic structure at the meta-GGA level.

% Additionally, for $c = 1$ the original Becke-Johnson potential is recovered \cite{Becke_JCP_2006}.

\subsection{The Heyd-Scuseria-Ernzerhof hybrid functional for solids}
\label{sec:ch_3_impl_HSE}

In Sec.~\ref{sec:ch_3_impl_LDA_TB-mBJ}, we described the determination of the XC functional using (i) the LDA, (ii) the PBE GGA, and (iii) the TB-mBJ meta-GGA. These approximations to the XC energy retain limitations and can produce systematic errors in the calculation of several key material properties \cite{Schimka_JCP_2011}. For example, the LDA (PBE) XC functional typically underestimates (overestimates) the lattice parameter by $1$ - $2$\%, and the error in PBE-calculated lattice parameters systematically increases with increasing ionic mass \cite{Schimka_JCP_2011}. To address limitations of this type, so-called generalised, or ``hybrid'', Kohn-Sham approaches have emerged \cite{Seidl_PRB_1996}. This generalised formalism involves introducing a portion of exact Hartree-Fock exchange into the XC energy. The HSE screened hybrid functional, introduced by Heyd, Scuseria and Ernzerhof \cite{Heyd_JCP_2003}, is one such approach. In general, hybrid XC functionals describe a number of chemical properties accurately but, typically, with the trade-off that the calculation of the Hartree-Fock exchange term is highly computationally expensive \cite{Heyd_JCP_2003}. The HSE functional is based on a screened Coulomb potential for the exchange interaction which partially mitigates this issue. In the HSE approach, the Coulomb kernel for the exchange energy is separated into short-range (SR) and long-range (LR) contributions following \cite{Caro_thesis_2013}

\begin{equation}
    \frac{1}{\vert \textbf{r} - \textbf{r}' \vert} = \underbrace{ \frac{\text{erfc} \left( \mu \, \vert \textbf{r} - \textbf{r}^{\prime} \vert \right)}{\vert \textbf{r} - \textbf{r}' \vert} }_{\text{SR}} + \underbrace{ \frac{\text{erf} \left( \mu \, \vert \textbf{r} - \textbf{r}^{\prime} \vert \right)}{\vert \textbf{r} - \textbf{r}' \vert} }_{\text{LR}} \, , \label{eq:ch_3_HSE_SR_LR_sep}
\end{equation}

where $\text{erf}$ and $\text{erfc}$ are the error and complementary error functions, respectively, and $\mu^{-1}$ is the distance at which the range separation is imposed. In the HSE approximation, only a portion of the SR exchange is Hartree-Fock (HF), the remainder of the exchange energy is determined via the PBE functional \cite{Caro_thesis_2013}. The HSE XC energy is expressed as \cite{Heyd_JCP_2003}

\begin{equation}
    E_{\scalebox{0.6}{\text{XC}}}^{\scalebox{0.6}{\text{HSE}}} = \alpha E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{HF,SR}}} (\mu) + (1-\alpha) E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{PBE,SR}}} (\mu) + E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{PBE,LR}}} (\mu) + E_{\scalebox{0.6}{\text{C}}}^{\scalebox{0.6}{\text{PBE}}} \, , \label{eq:ch_3_HSE_XC_energy}
\end{equation}

\noindent
where $\alpha$ is the exact-exchange mixing parameter determining how much exact-exchange is included in the calculation. We employ the typical value $\mu = 0.2$ \AA$^{-1}$ \cite{Caro_thesis_2013}. This approach yields an improved description of many chemical properties, particularly in comparison to the pure PBE functional \cite{Schimka_JCP_2011}. Notably, the HSE description of the XC energy effectively addresses the band gap underestimation problem that limits the applicability of the LDA, as well as many GGA and meta-GGA approaches (cf.~Sec.~\ref{sec:ch_3_impl_LDA_TB-mBJ}). However, the tendency of the PBE functional to overestimate lattice parameters and incorrectly describe atomisation energies is retained \cite{Schimka_JCP_2011}. To address these issues, the HSE functional for solids (HSEsol) was introduced by Schimka et al.~\cite{Schimka_JCP_2011}. The HSEsol XC functional takes the same form as the HSE functional in Eq.~\eqref{eq:ch_3_HSE_XC_energy}, but is based on the PBEsol functional for the semi-local exchange and correlation terms, where PBEsol is the revised PBE functional for solids \cite{Perdew_PRL_2008}. The use of the PBEsol XC functional provides an improved description of equilibrium properties of solids and surfaces \cite{Perdew_PRL_2008,Csonka_PRB_2009}. The use of the HSEsol hybrid functional has demonstrated an improved description of the chemical properties of solids vs.~other XC functionals, including the lattice parameters and enthalpies of formation \cite{Schimka_JCP_2011}, and notably the accurate description of band edge effective masses \cite{Laurien_PRB_2022}. Crucially, the HSEsol functional provides an improved description of these material properties while producing accurate band gaps \cite{Schimka_JCP_2011}. Like the BR mixing $c$ in the TB-mBJ exchange potential, the exact exchange mixing $\alpha$ in the HSE(sol) XC functional can be treated as an empirical parameter to fit to a target band gap. Therefore, in employing the HSEsol XC functional in our analysis of III-P band parameters in Chapter~\ref{ch:chapter6}, we apply this procedure and undertake a critical analysis of its performance vs.~experiment.

% $E_{\scalebox{0.6}{XC}}^{\scalebox{0.6}{\text{HSEsol}}}$

%\begin{equation}
%    E_{\scalebox{0.6}{XC}}^{\scalebox{0.6}{\text{HSEsol}}} = \alpha E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{HF,SR}}} (\omega) + (1-\alpha) E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{PBEsol,SR}}} %(\omega) + E_{\scalebox{0.6}{\text{X}}}^{\scalebox{0.6}{\text{PBEsol,LR}}} (\omega) + E_{\scalebox{0.6}{\text{C}}}^{\scalebox{0.6}{\text{PBEsol}}} \, . \label{eq:ch_3_HSEsol_XC_energy}
%\end{equation}

\subsection{Pseudopotentials and the projector augmented-wave method}
\label{sec:ch_3_impl_PAW}

% Pseudopotentials

A pseudopotential is an effective potential constructed to replace the strong Coulomb potential of the nucleus and the effects of the tightly bound core electrons by an effective ionic potential acting on the valence electrons beyond a certain cut-off distance from the nucleus, or ionic core \cite{Blochl_PRB_1994,Martin_book_2004}. This pseudopotential can then be substituted for the all-electron potential of Eq.~\eqref{eq:ch_3_KS_hamiltonian_martin}. In a psuedopotential calculation, the ``core'' consists of the nucleus and electrons occupying orbitals in filled shells, which are strongly localised around the nucleus \cite{Caro_thesis_2013}. Beyond the cut-off distance, outside the core region, the potential is that of the Coulomb interaction between the core and each of the valence electrons. Inside the core region the true potential, which is characterised by rapid spatial oscillations, is challenging to describe numerically, mandating very fine grid point spacing in real space calculations, or very high plane wave cut-off energy in reciprocal space calculations. Numerically, this true potential is replaced by a smooth function. Pseudopotentials are constructed so that the valence electron wave functions outside the core region are the same as those of an all-electron calculation. The use of pseudopotentials is strongly motivated by the fact that most chemical properties can be modelled accurately in calculations which only account for valence states. As outlined above, pseudopotentials additionally offer a numerically efficient platform to perform such calculations. The valence states typically have minimal interaction with the core states, a result of tightly bound core states possessing significantly lower energies. Notable exceptions to this are elements possessing ``shallow core'' or ``semi-core'' states, which possess energies large enough to influence chemical properties \cite{Caro_thesis_2013,Fiorentini_PRB_1993}.

% The PAW method

The PAW method, originally introduced by Bl{\"o}chl \cite{Blochl_PRB_1994}, has gained popularity for the description of psuedopotentials. The PAW method offers computational advantages compared to other approaches, such as norm-conserving and ultra-soft pseudopotentials \cite{Kresse_PRB_1999}. The PAW method considers a linear transformation $\widehat{\mathcal{T}}$ relating the all-electron valence wave function $\ket{\psi}$ and the smooth pseudo wave function $\ket{\psi_{\scalebox{0.6}{\text{PS}}}}$. $\widehat{\mathcal{T}}$ is equal to the identity operator $\widehat{I}$ outside of the core region, i.e.~$\ket{\psi}$ and $\ket{\psi_{\scalebox{0.6}{\text{PS}}}}$ are set equal outside of the core region, with $\widehat{\mathcal{T}} = \widehat{I} + \widehat{\mathcal{T}}_{R}$ defining the transformation operator $\widehat{\mathcal{T}_{R}}$ acting within the core region. $\widehat{\mathcal{T}}_{R}$ is specified by two sets of partial waves, the all-electron partial waves $\ket{\phi_{k}}$ and the pseudo partial waves $\ket{\phi_{k}^{\scalebox{0.6}{\text{PS}}}}$ which, within the core region, give

\begin{eqnarray}
    \ket{\psi} &=& \sum_{k} c_{k} \ket{\phi_{k}} \, , \label{eq:ch_3_all_electron_core} \\
%   
    \ket{\psi_{\scalebox{0.6}{\text{PS}}}} &=& \sum_{k} c_{k} \ket{\phi_{k}^{\scalebox{0.6}{\text{PS}}}} \, , \label{eq:ch_3_pseudo_core}
\end{eqnarray}

% $c_{k}^{\scalebox{0.6}{\text{PS}}} = \braket{ p_{k} \vert \phi_{k}^{\scalebox{0.6}{\text{PS}}} }$

\noindent
where the expansion coefficients $c_{k} = \braket{ p_{k} \vert \phi_{k}^{\scalebox{0.6}{\text{PS}}} }$ are determined by the projector $\bra{p_{k}}$, which resolves the identity via $\sum_{k} \ket{\phi_{k}^{\scalebox{0.6}{\text{PS}}}} \bra{p_{k}} = \widehat{I}$ within the core region. It follows that the transformation $\widehat{\mathcal{T}}$ can be written as

\begin{equation}
    \widehat{\mathcal{T}} = \widehat{I} + \sum_{k} \left( \ket{\phi_{k}} - \ket{\phi_{k}^{\scalebox{0.6}{\text{PS}}}} \right) \bra{p_{k}} \, . \label{eq:ch_3_transformation_paw}
\end{equation}

In the PAW method, the all-electron wave function $\ket{\psi}$ can then be obtained from the pseudo wave function $\ket{\psi^{\scalebox{0.6}{\text{PS}}}}$ as

\begin{equation}
    \ket{\psi} = \ket{\psi^{\scalebox{0.6}{\text{PS}}}} + \sum_{k} \left( \ket{\phi_{k}} - \ket{\phi_{k}^{\scalebox{0.6}{\text{PS}}}} \right) \braket{p_{k} \vert \psi^{\scalebox{0.6}{\text{PS}}}} \, . \label{eq:ch_3_all_electron_pseudo_wave_function_relation}
\end{equation}

Guidelines for the implementation of the all-electron operators as pseudo operators, and expressions for the charge density and total energy, are given by Bl{\"o}chl in Ref.~\cite{Blochl_PRB_1994}. Further technical details regarding the implementation of the PAW method are given by Kresse and Joubert in Ref.~\cite{Kresse_PRB_1999}.


\subsection{The plane wave formalism}
\label{sec:ch_3_impl_plane_waves}

As shown in Sec.~\ref{sec:ch_2_pwem}, formulating electronic structure calculations using plane waves is effective as a result of the periodic nature of crystalline materials. Expanding the Kohn-Sham wave functions in a plane wave basis is an approach employed by a number of DFT codes, including the VASP code employed in this work. The plane-wave expansion is expressed in three-dimensions as

\begin{equation}
     \psi_{n \textbf{k} } ( \textbf{r} ) = \frac{1}{ \sqrt{ \Omega } } \sum_{ \textbf{G} } a_{ n \textbf{k} \textbf{G} } \, e^{ i ( \textbf{k} + \textbf{G} ) \cdot \textbf{r} } \, ,
\end{equation}

\noindent
where $a_{n,\textbf{G}} (\textbf{k})$ are the Fourier expansion coefficients of the single-particle eigenstates, for which there is one eigenvalue $\epsilon_{n} (\textbf{k})$ per $\textbf{k}$-point per band $n$, we are operating in the momentum basis $\ket{\textbf{q}} = \ket{\textbf{k} + \textbf{G}}$, and note that $\langle \textbf{r} \vert \textbf{k} + \textbf{G} \rangle = \frac{ e^{ i ( \textbf{k} + \textbf{G} ) \cdot \textbf{r} } }{ \sqrt{ \Omega } }$, denoting the real-space representation of the plane wave basis states. In plane wave DFT calculations, the cut-off wave vector $G_{\text{cut}}$ determines the cut-off energy $E_{\text{cut}} = \frac{\hbar^{2}}{2m} G_{\text{cut}}^{2}$, which specifies the basis set via $\vert \textbf{G} + \textbf{k} \vert < G_{\text{cut}}$. As such, all plane waves with a kinetic energy smaller than $E_{\text{cut}}$ are included in the basis set. The choice of cut-off energy is analogous to the choice of the total number of plane waves in the plane wave expansion method (cf.~Sec.~\ref{sec:ch_2_pwem}). For any function $f(\textbf{r})$ sharing the same periodicity as the lattice, $f(\textbf{r}) = f(\textbf{r}+\textbf{R})$ (cf.~Sec.~\ref{sec:ch_2_k.p_theory_outline}), the Fourier transform $\widetilde{f} ( \textbf{G} )$ of $f(\textbf{r})$ is non-zero only when $\textbf{G}$ is a reciprocal lattice vector. The Fourier transform of the function $f(\textbf{r})$ is defined as

\begin{equation}
    \tilde{f} (\textbf{G}) = \frac{1}{\Omega} \int_{\Omega} f(\textbf{r}) e^{-i \textbf{G} \cdot \textbf{r}} \text{d}\textbf{r} \, , \label{eq:ch_3_fourier_transform}
\end{equation}

\noindent
where $\Omega = \textbf{a}_{1} \cdot ( \textbf{a}_{2} \times \textbf{a}_{3} )$ is the unit cell volume (cf.~Sec.~\ref{sec:ch_6_elastic_structural_lattice_param}). If the single-electron (Kohn-Sham) Hamiltonian $\widehat{H}$ is given by the sum of the kinetic energy operator and an effective potential $V_{\text{eff}} (\textbf{r})$, the Kohn-Sham equation can be written as \cite{Martin_book_2004}

\begin{equation}
    \sum_{\textbf{G}'} H_{\textbf{G},\textbf{G}'} (\textbf{k}) a_{n,\textbf{G}'} (\textbf{k}) = \epsilon_{n} (\textbf{k}) a_{n,\textbf{G}} (\textbf{k}) \, . \label{eq:ch_3_schrodinger_plane_wave}
\end{equation}

We note that the Fourier transform applied in obtaining Eq.~\eqref{eq:ch_3_schrodinger_plane_wave}, mirrors the plane wave expansion method introduced in the context of $\textbf{k} \cdot \textbf{p}$-based heterostructure calculations in Sec.~\ref{sec:ch_2_pwem}. These approaches to solving the relevant Schr\"{o}dinger equations -- i.e.~the coupled multi-band envelope function Schr\"{o}dinger equations in $\textbf{k} \cdot \textbf{p}$ calculations, or the Kohn-Sham equations in DFT -- are formally equivalent, and lead to a reciprocal space representation of Schr{\"o}dinger's equation. The matrix elements $H_{\textbf{G},\textbf{G}'}$ of the Kohn-Sham Hamiltonian are given by

\begin{equation}
    H_{\textbf{G},\textbf{G}'} (\textbf{k}) = \frac{\hbar^{2}}{2m_{0}} \vert \textbf{k} + \textbf{G} \vert^{2} \delta_{\textbf{G},\textbf{G}'} + \widetilde{V}_{\text{eff}} (\textbf{G} - \textbf{G}') \, , \label{eq:ch_3_eff_hamiltonian_matrix_elements}
\end{equation}

\noindent
where $\widetilde{V}_{\text{eff}} (\textbf{G} - \textbf{G}')$ is the Fourier transform of the effective potential (cf.~Sec.~\ref{sec:ch_3_kohn_sham}). Convergence in plane wave DFT calculations is determined, in part, by the size of the corresponding basis set \cite{Caro_thesis_2013}. Computing the electronic contribution to the total energy mandates including all contributions from filled (valence) bands at all wave vectors $\textbf{k}$ in the Brillouin zone. The Kohn-Sham equations must be solved at wave vectors \textbf{k} throughout the Brillouin zone and an integral with respect to \textbf{k} must then be performed. Therefore, as Eq.~\eqref{eq:ch_3_schrodinger_plane_wave} must be solved at each $\textbf{k}$-point, the evaluation of any function accurately becomes computationally expensive in realistic calculations with large numbers of $\textbf{k}$-points. Given a \textbf{k}-dependent function $\eta_{n} (\textbf{k})$, the integral, where the summation takes place over energy bands $n$, is

\begin{eqnarray}
    \int \eta_{n} \, \text{d}\textbf{k} &=& \frac{1}{(2\pi)^{3}} \sum_{n} \int_{\text{BZ}} \eta_{n} (\textbf{k}) f_{n,\textbf{k}} \text{d}\textbf{k} \, , \nonumber \\
    &\rightarrow& \frac{1}{ \left( 2 \pi \right)^{3} } \frac{1}{\Omega} \sum_{ \textbf{k} } w_{\textbf{k}} \, \eta_{n}(\textbf{k}) \, f_{n,\textbf{k}} \, , \label{eq:ch_3_parameter_plane_wave}
\end{eqnarray}

%\begin{equation}
%    \int_{\text{BZ}} \frac{ \text{d} \textbf{k} }{ ( 2 \pi )^{3} } \to \frac{1}{ \Omega } \sum_{ \textbf{k} } w_{\textbf{k}}
%\end{equation}

\noindent
where the occupancy $f_{n,\textbf{k}}$, in the zero-temperature limit, is given by 1 (occupied) or 0 (unoccupied), and $w_{\textbf{k}}$ is a weighting factor related to the distribution of $\textbf{k}$ states chosen. In the finite-temperature case, the occupancies can vary continuously between 0 and 1. The accuracy of calculating any function following Eq.~\eqref{eq:ch_3_parameter_plane_wave} is determined by the density and choice of the weighted sum over $\textbf{k}$-points.

As the number of $\textbf{k}$-points included in a calculation increases, the accuracy of the numerical approximation to the integral in Eq.~\eqref{eq:ch_3_parameter_plane_wave} will improve \cite{Caro_thesis_2013}. However, calculations using large numbers of $\textbf{k}$-points quickly become computationally prohibitive, meaning that there must be careful consideration of the $\textbf{k}$-points employed in DFT calculations. For energies near the Fermi level, the occupancies will, in general, vary for the same band at different $\textbf{k}$-points. For narrow-gap semiconductors and metals, a large number of $\textbf{k}$-points may be required to reconstruct the Fermi surface \cite{Fu_PRB_1983}. This is not an issue for the relatively wide gap semiconductors considered in this Thesis, but can pose computational challenges for convergence vs.~the number of \textbf{k}-points employed in the BZ integration in other cases. The smearing method, originally introduced by Fu and Ho \cite{Fu_PRB_1983}, is one technique to accelerate this convergence with $\textbf{k}$-points. This approach ``smears'' the energies of each of the eigenstates using, e.g., a Gaussian distribution. The Fermi energy is then calculated from the density of states given by the smeared energies. Then, the partial occupancy of each band equals the portion of its smeared value lying below the calculated Fermi energy \cite{Fu_PRB_1983}. This mandates the total energy $E$ being replaced by a generalised free energy $F$, that latter containing an extra term that depends on the partial occupancies \cite{Kresse_PRB_1996}, resembling the entropy of the finite-temperature case \cite{Caro_thesis_2013}.

% As a result, in our DFT calculations we select specially constructed $\textbf{k}$-points grids, with careful consideration in selecting the density of these $\textbf{k}$-point grids and the cut-off wave vector (energy) $G_{\text{cut}}$ ($E_{\text{cut}}$).

% Additionally, Eq.~\eqref{eq:ch_3_parameter_average_plane_wave} employs an integral over the BZ instead of a summation as for macroscopic crystals the allowed values of $\textbf{k}$ are quasi-continuous \cite{Ibach_book_2003,Martin_book_2004}.

We note that, in numerical Kohn-Sham DFT calculations, the Schr{\"o}dinger equation is solved self-consistently via an iterative procedure. Beginning from an initial charge density and iteratively optimising towards convergence, the Kohn-Sham Hamiltonian is specified and the single-particle eigenstates computed, with these eigenstates then employed to construct an updated charge density (cf.~Eq.~\eqref{eq:ch_3_KS_charge_density}). This updated charge density then serves as the input for the next iteration, and the procedure is iterated until the difference in total energy between two consecutive iterations has converged to a desired precision -- i.e.~achieving so-called charge self-consistency. Upon the completion of this self-consistent cycle the final single-particle eigenstates (Kohn-Sham orbitals) define a converged charge density which, due to the Hohenberg-Kohn theorems, uniquely determines the electronic ground state \cite{Hohenberg_PR_1964}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Section 3: Conclusions %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Conclusions}
\label{sec:ch_3_conclusions}

In this Chapter we have presented the theoretical methods employed in our analysis of the electronic band parameters of the III-P compounds AlP, GaP and InP in Chapter~\ref{ch:chapter6}. We began in Sec.~\ref{sec:ch_3_kohn_sham} by presenting an overview of the Kohn-Sham formalism of DFT. We described the approach to compute the total energy of a many-electron system in this formalism, where having used the Born-Oppenheimer approximation to freeze the nuclei and separate the electrons and nuclei, and their interactions, the many-electron problem is tackled using the single-electron Kohn-Sham approach. We described how, under the application of appropriate approximations, the accuracy of the implementation of the Kohn-Sham formalism can be largely determined by the XC energy. The XC energy is defined as the energy difference between the exact many-body interacting system and the non-interacting electron density, described by a set of independent wave functions. In Secs.~\ref{sec:ch_3_impl_LDA_TB-mBJ} and \ref{sec:ch_3_impl_HSE}, we then provided detail on the treatment of the XC energy as a functional of the electron density. In Sec.~\ref{sec:ch_3_impl_LDA_TB-mBJ}, we began with the description of the LDA, an XC functional corresponding to a UEG where the electron density is treated as constant. We then described the GGA and meta-GGA, which go beyond the LDA by including semi-local effects via spatial derivatives of the electron density. We described the TB-mBJ exchange potential in which the Coulomb potential created by the exchange hole is modelled via the BR potential. The TB-mBJ scheme is a potential only approach and must therefore be employed in conjunction with the LDA correlation, but delivers improved description of the electronic structure vs.~LDA alone. Notably the TB-mBJ approach overcomes the band gap underestimation problem present in the LDA. In Sec.~\ref{sec:ch_3_impl_HSE}, we outlined the hybrid HSEsol XC functional, a hybrid XC functional that combines the PBEsol functional with a fraction $\alpha$ of the exact Hartree-Fock exchange energy. In Sec.~\ref{sec:ch_3_impl_HSE}, we described that the HSEsol functional provides an accurate description of the chemical properties of solids in Kohn-Sham formalised DFT, including providing an accurate description of the band edge effective masses and the band gaps of semiconductors and insulators. Next, in Sec.~\ref{sec:ch_3_impl_PAW} we described pseudopotentials, and the PAW method, which is a robust and efficient approach to compute the all-electron wave function, that represents an extension of the commonly employed pseudopotential approach. In Sec.~\ref{sec:ch_3_impl_plane_waves}, we then outline the plane wave formalism for DFT calculations, and how this technique offers a computationally efficient basis for DFT calculations, akin to benefits offered by the plane wave expansion method seen in Chapter~\ref{ch:chapter2}. The theoretical methods and technical details outlined in this Chapter will underpin our analysis of III-P band parameters in Chapter~\ref{ch:chapter6}.