%****************************************************************
\section{Schr\"odinger equation}
%****************************************************************
In this memoir, we consider atomic and molecular quantum systems (i.e.~systems composed by nuclei and electrons) within the Born-Oppenheimer approximation \cite{Szabo, Helgaker}.
This means that we neglect the kinetic energy of the nuclei and treat the nuclear coordinates as parameters.
We therefore concentrate our attention on the electronic degrees of freedom.
Unless otherwise stated, atomic units are used throughout this memoir.

A chemical system is completely defined at a time $t$ by its electronic wave function $\Psi(\bX,t)$, solution of the time-dependent Schr\"odinger equation
\begin{equation} 
\label{eq:time-dependent-schrodinger-equation}
	i \pdv{\Psi(\bX,t)}{t} = \HOp \Psi(\bX,t),
\end{equation}
where $\HOp$ is the so-called Hamiltonian and $\bX = (\bx_1,\dots,\bx_n) = (\bs,\bR)$ is a composite coordinate vector gathering the spin coordinates $\bs = (s_1,\ldots,s_n)$ and spatial coordinates $\bR = (\br_1,\ldots,\br_n)$ of the $n$ electrons.
 
In the case of a stationary system, the time-independent Schr\"odinger equation reads
\begin{equation}
\label{eq:time-independent-schrodinger-equation}
	\HOp \Psi(\bX,t) = E\,\Psi(\bX,t),
\end{equation}
where $E$ is the energy of the system and the non-relativistic Hamiltonian is explicitly given by
\begin{equation} 
\label{eq:Hamiltonian}
\begin{split} 
	\HOp	& = \TeOp + \VenOp + \VeeOp + \VnnOp
			\\
			& = - \sum_i^n \frac{\nabla^2_i}{2}
			- \sum_i^n \sum_A^\nuc \frac{Z_A}{\abs{\brA - \bri}}
			+ \sum_{i < j}^n \frac{1}{\abs{\bri - \brj}} 
			+ \sum_{A<B}^\nuc \frac{Z_A Z_B}{\abs{\brB - \brA}},
\end{split}
\end{equation}
where $\nabla_i^2$ is the Laplace operator associated with the $i$th electron, $\brA$ and $Z_A$ are the nuclear coordinates and charge of nucleus $A$.
The first term $\TeOp$ is the kinetic energy operator of the electrons, the next term $\VenOp$ corresponds to the Coulombic attraction between electrons and nuclei, while the last two terms corresponds to the interelectronic ($\VeeOp$) and internuclear ($\VnnOp$) Coulombic repulsions, respectively.
Note that the last term $\VnnOp$ is a constant as, within the Born-Oppenheimer approximation, it does not depend on the electronic coordinates, and will be omitted for the sake of clarity.

%****************************************************************
\section{Hartree-Fock approximation}
%****************************************************************
Within the Hartree-Fock (HF) approximation \cite{Szabo}, the electronic wave function $\PsiHF$ is written as a Slater determinant of $n$ spin orbitals
\begin{equation} 
\label{eq:Psi-HF}
	\PsiHF(\bx_1,\ldots,\bx_n)  = \frac{1}{\sqrt{n!}} 
	\begin{vmatrix}
		\sMO{1}(\bx_1) 	& \sMO{2}(\bx_1)	& \ldots 	& \sMO{n}(\bx_1)	\\
		\sMO{1}(\bx_2) 	& \sMO{2}(\bx_2) 	& \ldots 	& \sMO{n}(\bx_2) 	\\
		\vdots 		& \vdots		& \ddots	& \vdots 		\\ 
		\sMO{1}(\bx_n) 	& \sMO{2}(\bx_n) 	& \ldots 	& \sMO{n}(\bx_n) 	\\
	\end{vmatrix}.
\end{equation}
Each spin orbital $\sMO{i}(\bx)$ is a product of a spin part $\omega(s)$ and a spatial part $\MO{i}(\br)$ --- also known as a molecular orbital (MO) ---
\begin{equation}
	\sMO{i}(\bx) = \omega(s) \MO{i}(\br),
\end{equation}
where
\begin{equation}
	\omega(s) = 
	\begin{cases}
		\alpha(s),	&	\text{for spin-up electrons,}	
		\\
		\beta(s),	&	\text{for spin-down electrons.}
	\end{cases}
\end{equation}
The spin orbitals form an orthonormal set, i.e.
\begin{equation} 
	\braket{\sMO{i}}{\sMO{j}} = \delta_{ij},
\end{equation}
where 
\begin{equation} 
	\delta_{ij} = 
	\begin{cases}
		1, & \text{if $i = j$},\\
		0, & \text{otherwise},
	\end{cases}
\end{equation}
is the Kronecker delta \cite{NISTbook}. 

The HF energy is defined as
\begin{equation} 
	\EHF = \mel{\PsiHF}{\HOp}{\PsiHF},
\end{equation}
and yields the following expression:
\begin{equation}
\begin{split}
	\EHF 	& = \sum_i^n \mel{\sMO{i}(\br_1)}{\HcOp}{\sMO{i}(\br_1)}
		\\
		& + \sum_{i<j}^n \qty[ 
		\mel{\sMO{i}(\br_1)\sMO{j}(\br_2)}{\ree^{-1}}{\sMO{i}(\br_1)\sMO{j}(\br_2)}
		- \mel{\sMO{i}(\br_1)\sMO{j}(\br_2)}{\ree^{-1}}{\sMO{j}(\br_1)\sMO{i}(\br_2)} ],
\end{split}
\end{equation}
where the so-called core Hamiltonian (i.e.~the one-electron part of the electronic Hamiltonian) is defined as
\begin{equation}
	\HcOp = \TeOp + \VenOp.
\end{equation}
We define the Fock operator as
\begin{gather}
	\FOp \sMO{i}(\br_1) = \MOev{i} \sMO{i}(\br_1),
	\\
	\FOp(\br_1) = \HcOp(\br_1) + \sum_i^n \qty[ \JOp_i(\br_1) - \KOp_i(\br_1) ],
\end{gather}
where $\JOp_i(\br_1)$ and $\KOp_i(\br_1)$ are the Coulomb and exchange operators respectively:
\begin{subequations}
\begin{align} 
	\JOp_i(\br_1) \sMO{j}(\br_1) & = \sMO{j}(\br_1) \int \sMO{i}(\br_2) \ree^{-1} \sMO{i}(\br_2) d\br_2,
	\\
	\KOp_i(\br_1) \sMO{j}(\br_1) & = \sMO{i}(\br_1) \int \sMO{i}(\br_2) \ree^{-1} \sMO{j}(\br_2) d\br_2.
\end{align}
\end{subequations}
In the following, we adopt the restricted HF (RHF) formalism which means that we assume that the spatial part of the spin orbital is independent of the spin state of the electron occupying this orbital \cite{Szabo}.
Moreover, unless otherwise stated, the systems treated here are closed-shell systems, i.e.~each MO is doubly occupied by one spin-up and one spin-down electron.

%****************************************************************
\subsection{Roothaan-Hall equations}
%****************************************************************
Within the LCAO approximation, we expand each MO as a linear combination of $N$ atomic orbitals (AOs), such as
\begin{equation} 
\label{eq:LCAO}
	\MO{i}(\br) = \sum_\mu \cMO{\mu}{i} \AO{\mu}(\br).
\end{equation}
In practice, the AOs $\AO{\mu}(\br)$ are usually chosen as cartesian Gaussian functions due to their computational convenience.
However, other choices (such a Slater functions) are possible depending on the type of systems and the target accuracy.
We will come back to this particular point later in this memoir.

In the AO basis, we have 
\begin{align} 
\label{eq:FAO}
	\FkEl{\mu}{\nu} = \mel{\AO{\mu}}{\FOp}{\AO{\nu}}
			\equiv \mel{\mu}{\FOp}{\nu}
			& = \HcEl{\mu}{\nu} + \sum_{\lambda \sigma} \PEl{\lambda}{\sigma} 
			\qty[ \braket{\mu \lambda}{\nu \sigma} - \frac{1}{2} \braket{\mu\lambda}{\sigma \nu} ],
\end{align}
with 
\begin{gather}
	\HcEl{\mu}{\nu} = \mel{\mu}{\HcOp}{\nu},
	\\
	\braket{\mu \lambda}{\nu \sigma} = \iint \AO{\mu}(\br_1) \AO{\lambda}(\br_2) \ree^{-1} \AO{\nu}(\br_1) \AO{\sigma}(\br_2) d\br_1 d\br_2,
\end{gather}
and where the density matrix is defined as
\begin{equation} 
	\PEl{\mu}{\nu} = 2 \sum_i^\occ \cMO{\mu}{i} \cMO{\nu}{i}.
\end{equation}
The HF electronic energy of the system is then given by
\begin{equation} 
\label{eq:E-elec}
	\EHF = \sum_{\mu \nu} \PEl{\mu}{\nu} \HcEl{\mu}{\nu}
	+ \frac{1}{2} \sum_{\mu \nu \lambda \sigma}^N \PEl{\mu}{\nu} \PEl{\lambda}{\sigma} G_{\mu \nu \lambda \sigma}.
\end{equation}
In matrix form, the Fock matrix $\FkMat$ can be decomposed as
\begin{equation}
	\FkMat = \HcMat + \GMat,
\end{equation}
where
\begin{align}
	\GEl{\mu}{\nu} 				& = \sum_{\lambda \sigma} \PEl{\lambda}{\sigma} \GEl{\mu\nu}{\lambda\sigma},
	&
	\GEl{\mu\nu}{\lambda\sigma} & = \braket{\mu \lambda}{\nu \sigma} - \frac{1}{2} \braket{\mu\lambda}{\sigma \nu}.
\end{align}
The stationnarity of the energy with respect to the coefficients $\cMO{\mu}{i}$ yields the Roothaan-Hall equations:
\begin{equation}
\label{eq:Roothaan-Hall}
	\sum_\nu \FEl{\mu}{\nu} \cMO{\nu}{i} = \sum_\nu \SEl{\mu}{\nu} \cMO{\nu}{i} \MOev{i},
\end{equation}
or in matrix form
\begin{equation} 
\label{eq:Roothaan-Hall-mat}
	\FkMat \CMat = \SMat \CMat \MOevMat,
\end{equation}
where the elements of the overlap matrix $\SMat$ are given by
\begin{equation}
	\SEl{\mu}{\nu} = \braket{\mu}{\nu}.
\end{equation}
The coefficient matrix $\bC$ gathers the MO coefficients $\cMO{\mu}{i}$, while the diagonal matrix $\MOevMat$ gathers the MO energies $\MOev{i}$.
We introduce the orthogonalization matrix $\bX$ such as
\begin{equation} 
	 \XMat^\dag \SMat \XMat = \IdMat
\end{equation}
in order to work in an orthogonal AO basis (where $\IdMat$ is the identity matrix).
There are two main orthogonalisation methods, namely the L\"owdin orthogonalisation for which $\XMat = \SMat^{-1/2}$ and the canonical orthogonalisation for which $\XMat = \UMat \sMat^{-1/2}$ (where $\UMat$ and $\sMat$ are the eigenvectors of eigenvalues matrices of $\SMat$, respectively).
Nowadays, the usual procedure consists in performing a singular value decomposition (SVD) of the overlap matrix $\SMat$. 
This procedure is efficient, numerically stable and allows to remove the linear dependencies which might be present in the AO basis.

Rotating the Fock matrix $\FMat$ into the orthogonal basis yields
\begin{equation}
\label{eq:Fprime}
	\FkMat^{\prime} \CMat^{\prime} = \CMat^{\prime} \MOevMat,
\end{equation}
where
\begin{equation}
	\FkMat^{\prime} = \XMat^{\dag} \FMat \XMat.
\end{equation}
The matrices $\CMat^{\prime}$ and $\MOevMat$ can be determined by a straightforward diagonalisation of Eq.~\eqref{eq:Fprime}, and the matrix $\CMat$ is obtained by back-transforming the eigenvectors in the original basis: 
\begin{equation}
	\CMat = \XMat \CMat^{\prime}.
\end{equation}

%****************************************************************
\subsection{Self-consistent field calculation}
%****************************************************************
In order to obtain the MO coefficients $\CMat$, one must diagonalise the Fock matrix $\FkMat$.
However, this matrix does depend on the MO coefficients itself.
Therefore, one must employ an iterative procedure called self-consistent field (SCF) method.
The SCF algorithm is described below:
\begin{enumerate}
	\item Obtain an estimate of the density matrix $\PMat$.
	\item Build the Fock matrix: $\FMat = \HcMat + \GMat$.
	\item Transform the Fock matrix in the orthogonal matrix: $\FMat^{\prime} = \XMat^{\dag} \FMat \XMat$.
	\item Diagonalize $\FMat^{\prime}$ to obtain $\CMat^{\prime}$ and $\MOevMat$.
	\item Back-transform the MOs in the original basis: $\CMat = \XMat \CMat^{\prime}$.
	\item Compute the new density matrix $\PMat = \CMat \CMat^{\dag}$, as well as the HF energy: 
	\begin{equation}
		\EHF = \frac{1}{2} \Tr{\PMat \qty( \HcMat + \FkMat )}.
	\end{equation}
	\item Convergence test. If not satisfied, go back to $2.$
\end{enumerate}

Unfortunately, the HF method cannot be used to obtain the exact energy of the system even in the complete basis set (CBS) limit due to the approximate treatment of the electron-electron interaction. 
Within the HF method, this interaction is averaged over all the electrons.
In other word, a given electron ``feels'' the averaged repulsion of the $n-1$ remaining electrons (mean-field approach).
We will see in the next section how one can go beyond the HF approximation. 

%****************************************************************
\section{Post Hartree-Fock methods}
%****************************************************************
The correlation energy $\Ec$ is defined as the error in the HF approximation, i.e.~the energy difference between the exact energy and the energy calculated within the HF approximation:
\begin{equation} 
\label{eq:Ec}
	\Ec = E - \EHF.
\end{equation}
Thanks to the variational principle, $\Ec$ is always a negative quantity.
The purpose of post HF methods is to recover some or all of the correlation energy \cite{JensenBook, CramerBook}.

%****************************************************************
\subsection{Configuration interaction methods}
%****************************************************************
One of the most conceptually simple (albeit expensive) approach to recover a large fraction of the correlation energy is the configuration interaction (CI) method.
The general idea is to expand the wave function as a linear combination of ``excited'' determinants.
These excited determinants are built by promoting electrons from occupied to unoccupied (virtual) MOs usually based on the HF orbitals, i.e.
\begin{equation} 
\label{eq:Psi-CI}
	\PsiCI	= 					\cCI{0}{} 	\PsiHF 
		+ \sum_i^\occ \sum_a^\virt 		\cCI{i}{a} 	\ExDet{i}{a} 
		+ \sum_{ij}^\occ \sum_{ab}^\virt 	\cCI{ij}{ab} 	\ExDet{ij}{ab}
		+ \sum_{ijk}^\occ \sum_{abc}^\virt 	\cCI{ijk}{abc} 	\ExDet{ijk}{abc}
		+ \ldots,
\end{equation}
where $\ExDet{i}{a}$, $\ExDet{ij}{ab}$ and $\ExDet{ijk}{abc}$ are singly-, doubly- and triply-excited determinants.
$\ExDet{ij}{ab}$ corresponds to the excitations of two electrons from the occupied spinorbitals $i$ and $j$ to virtual spinorbitals $a$ and $b$. 
It is easy to show that the CI energy
\begin{equation}
	\ECI = \mel{\PsiCI}{\HOp}{\PsiCI}
\end{equation}
is an upper bound to the exact energy of the system.

When all possible excitations are taken into account, the method is called full CI (FCI) and it recovers the entire correlation energy for a given basis set.
Albeit elegant, FCI is very expensive due to the exponential increase of the number of excited determinants.
For example, when only singles and doubles are taken into account, the method is called CISD. 
It recovers an important chunk of the correlation.
However, it has the disadvantage to be size-inconsistent.

A method is said to be size-consistent if the correlation energy of two non-interaction systems is egal to twice the correlation energy of the isolated system.
Size-extensivity means that the correlation energy grows linearly with the system size.

%****************************************************************
\subsection{Density-functional theory}
%****************************************************************
Density-functional theory (DFT) is based on two theorems known as the Hohenberg-Kohn (HK) theorems \cite{Hohenberg64}, which states that it exists a non-interacting reference system with an electronic density $\rho(\br)$ equal to the real, interaction system. 
The first theorem proves the existence of a one-to-one mapping between the electron density and the external potential, while the second HK theorem guarantees the existence of a variational principle for the ground-state electron density.

%****************************************************************
\subsection{Kohn-Sham equations}
%****************************************************************
Present-day DFT calculations are almost exclusively done within the so-called Kohn-Sham (KS) formalism, which corresponds to an exact dressed one-electron theory \cite{Kohn65}.
In analogy to the HF theory, the electrons are treated as independent particles moving in the average field of all others but now with exchange an correlation included by virtue of an ``exchange-correlation'' functional.

Following the work of Kohn and Sham \cite{Kohn65}, we introduce KS orbitals $\sMO{i}(\br)$, and the energy can be decomposed as 
\begin{equation}
	\EKS \qty[\rho(\br)] 	= \Ts \qty[\rho(\br)] 
				+ \Ene \qty[\rho(\br)]
				+ J \qty[\rho(\br)]
				+ \Exc \qty[\rho(\br)],
\end{equation}
where 
\begin{equation}
	\Ts \qty[\rho(\br)] = - \frac{1}{2} \sum_i^\occ \mel{\sMO{i}}{\nabla_i^2}{\sMO{i}}
\end{equation}
is the non-interacting kinetic energy, 
\begin{equation}
	\Ene \qty[\rho(\br)] = - \sum_A^\nuc \int \frac{Z_A \rho(\br)}{\abs{\brA - \br}} d\br
\end{equation}
is the electron-nucleus attraction energy,
\begin{equation}
	J \qty[\rho(\br)] = \frac{1}{2} \iint \frac{\rho(\br_1) \rho(\br_2)}{\abs{\br_1 - \br_2}} d\br_1 d\br_2
\end{equation}
is the classical electronic repulsion, and the one-electron density is
\begin{equation}
        \rho(\br) = \sum_i^\occ \abs{\sMO{i}(\br)}^2.
\end{equation}
The exchange-correlation energy 
\begin{equation}
        \Exc \qty[\rho(\br)] = \qty{ T \qty[\rho(\br)] - \Ts \qty[\rho(\br)] } + \qty{ \Eee \qty[\rho(\br)] -  J \qty[\rho(\br)]}
\end{equation}
is the sum of two terms: one coming from the difference between the exact kinetic energy $T\qty[\rho(\br)]$ and the non-interacting kinetic energy $\Ts \qty[\rho(\br)]$, and the other one coming from the difference between the exact interelectronic repulsion $\Eee \qty[\rho(\br)]$ and the classical Coulomb repulsion $J \qty[\rho(\br)]$.
Here, we will only consider the second term as the ``kinetic'' correlation energy is usually much smaller than its ``Coulomb'' counterpart.
Also, as it is usually done, we will split the exchange-correlation energy as a sum of an exchange and correlation components, i.e.
\begin{equation}
        \Exc \qty[\rho(\br)] = \Ex \qty[\rho(\br)] + \Ec \qty[\rho(\br)].
\end{equation}

Similarly to the Roothaan-Hall equations, the condition of stationarity of the KS energy with respect to the electron density
\begin{equation}
	\frac{\delta E \left[ \rho (\mathbf{r}) \right]}{\delta \rho (\mathbf{r})} = \mu
\end{equation}
(where $\mu$ is the chemical potential) yields the KS equations
\begin{equation}
	 \qty[ - \frac{\nabla_{\br}^2}{2}
	- \sum_A^\nuc  \frac{Z_A}{\abs{\br - \brA}} 
	+ \int \frac{\rho(\br^{\prime})}{\abs{\br - \br^{\prime}}} d\br^{\prime} 
	+ \fdv{\Exc\qty[\rho(\br)]}{\rho(\br)} ] \sMO{i}(\br) 
	= \MOev{i} \sMO{i}(\br),
\end{equation}
which can be re-written as
\begin{equation}
	\FOp_\text{KS} \sMO{i}(\br) = \MOev{i} \sMO{i}(\br).
\end{equation}
These equations are solved iteratively, just like the HF equations, by expanding the KS MOs in a AO basis, yielding
\begin{equation}
	\FkMat_\text{KS}  \CMat = \SMat \CMat \MOevMat.
\end{equation}

%****************************************************************
\subsection{Exchange-correlation functionals}
%****************************************************************
Due to its moderate computational cost and its reasonable accuracy, KS DFT \cite{Hohenberg64, Kohn65} has become the workhorse of electronic structure calculations for atoms, molecules and solids \cite{ParrYang}.
To obtain accurate results within DFT, one only requires the exchange and correlation functionals, which can be classified in various families depending on their physical input quantities \cite{Becke14, Yu16}.
These various types of functionals are classified by the Jacob's ladder of DFT \cite{Perdew01, Pewdew05} (see Fig.~\ref{fig:Jacob-ladder}).

%%% FIGURE 1 %%%
\begin{figure}
	\centering
	\includegraphics[width=0.25\linewidth]{../Chapter1/fig/fig1}
\caption{
\label{fig:Jacob-ladder}
Jacob's ladder of DFT. 
$\rho$, $x$, $\tau$ and $\Ex^\text{HF}$ are the electron density, the reduced gradient, the kinetic energy density, and the HF exchange energy, respectively.}
\end{figure}
%%%

\begin{itemize}
	\item The local-density approximation (LDA) sits on the first rung of the Jacob's ladder and only uses as input the electron density $\rho$.
	The oldest and probably most famous LDA functional is the Dirac exchange functional (D30) \cite{Dirac30} based on the uniform electron gas (UEG) \cite{WIREs16}.
        Based on the work of Ceperley and Alder who used quantum Monte Carlo calculations (see below) to determine the correlation energy of the UEG with respect to the density \cite{Ceperley80},
	Vosko, Wilk and Nusair (VWN) proposed a LDA correlation functional by fitting their data \cite{VWN80}.

	\item The generalized-gradient approximation (GGA) corresponds to the second rung and adds the gradient of the electron density $\nabla \rho$ as an extra ingredient.
	The well-known B88, G96, PW91 and PBE exchange functionals are examples of GGA exchange functionals \cite{B88, G96, PW91, PBE}.
	Probably the most famous GGA correlation functional is LYP \cite{LYP}, which gave birth to the GGA exchange-correlation functional BLYP \cite{Becke88b} by combination with B88.

	\item The third rung is composed by the so-called meta-GGA (MGGA) functionals \cite{Sala16} which uses, in addition to $\rho$ and $\nabla \rho$, the kinetic energy density 
	\begin{equation}
		\tau = \sum_i^\occ \abs{\nabla\sMO{i} }^2.
	\end{equation}
	The M06-L functional from Zhao and Truhlar \cite{M06L}, the mBEEF functional from Wellendorff et al.~\cite{mBEEF} and the SCAN \cite{SCAN} and MS \cite{MS0,MS1_MS2} family of functionals from Sun et al.~are examples of widely-used MGGA functionals.

	\item The fourth rung (hyper-GGAs or HGGAs) includes the widely-used hybrid functionals, introduced by Becke in 1993 \cite{Becke93}, which add a certain percentage of HF exchange.
	Example of such functionals are B3LYP \cite{Becke93}, B3PW91 \cite{PW92, Becke93, Perdew96}, BH\&HLYP \cite{Becke93b} or PBE0 \cite{PBE0}.
	Hybrids functionals are known for their accuracy in electronic structure theory.
	However, they are more computationally expensive than LDA or GGA functionals due to the calculation of the costly HF exchange.
	
	\item The fifth rung includes double hybrids and RPA-like functionals but we will not be discussing such types of functionals in the present memoir.
\end{itemize}

%****************************************************************
\section{Quantum Monte Carlo methods}
%****************************************************************
%-----------------------------------------------
\subsection{Variational Monte Carlo}
%-----------------------------------------------
In the VMC method, the expectation value of the Hamiltonian with respect to a trial wave function is obtained using a stochastic integration technique. 
Within this approach a variational trial wave function $\PsiT(\bR)$ is introduced, and one then calculates its variational energy
\begin{equation}
	\EVMC = \frac{\int \PsiT(\bR) \HOp \PsiT(\bR) d\bR}{\int \PsiT(\bR)^2 d\bR},
\end{equation}
using the Metropolis Monte Carlo method of integration \cite{Umrigar99}.
The resulting VMC energy is an upper bound to the exact ground-state energy, within the statistical Monte Carlo error.  
Unfortunately, any resulting observables are biased by the form of the trial wave function, and the method is therefore only as good as the chosen $\PsiT$.

%--------------------------------------------
\subsection{Diffusion Monte Carlo}
%--------------------------------------------
DMC is a stochastic projector technique for solving the many-body Schr\"odinger equation \cite{Kalos74, Ceperley79, Reynolds82}.
Its starting point is the time-dependent Schr\"odinger equation in imaginary time
\begin{equation} 
\label{eq:DMC}
	-\pdv{\Psi(\bR,\tau)}{\tau} = (\HOp - S) \Psi(\bR,\tau).
\end{equation}
For $\tau \to \infty$, the steady-state solution of Eq.~\eqref{eq:DMC} for $S$ close to the ground-state energy is the ground-state $\Psi(\bR)$ \cite{Kolorenc11}.
DMC generates configurations distributed according to the product of the trial and exact ground-state wave functions.  
If the trial wave function has the correct nodes, the DMC method yields the exact energy, within a statistical error that can be made arbitrarily small by increasing the number of Monte Carlo steps. 
Thus, as in VMC, a high quality trial wave function is essential in order to achieve high accuracy \cite{Umrigar93, Huang97}.

%------------------------------------------
\subsection{Trial wave functions}
%------------------------------------------
Within QMC, trial wave functions are usually defined as \cite{Huang97, Drummond04, LopezRios12}
\begin{equation}
\label{eq:Psi-trial}
 	\PsiT(\bR) = e^{\Js(\bR)} \sum_I \cCI{I}{} \sdet_I^{\uparrow}(\br^{\uparrow}) \, \sdet_I^{\downarrow}(\br^{\downarrow}),
\end{equation}
where $D_I^{\sigma}$ are determinants of the spin-$\sigma$ electrons.
The fermionic nature of the wave function is imposed by a single- or multi-determinant expansion of Slater determinants made of HF or KS MOs.
$\Js(\bR)$ is called the Jastrow factor and $e^{\Js(\bR)}$ is a nodeless function.
Hence, the nodes of $\PsiT$ are completely determined by the determinantal part of the trial wave function.

%---------------------------------------------------
\subsection{Fixed-node approximation}
%---------------------------------------------------
Considering an antisymmetric (real) electronic wave function $\Psi(\bR)$, the nodal hypersurface (or simply ``nodes'') is a $(n\, D-1)$-dimensional manifold defined by the set of configuration points $\bN$ for which $\Psi(\bN)=0$.
The nodes divide the configuration space into nodal cells or domains which are either positive or negative depending on the sign of the electronic wave function in each of these domains.
In recent years, strong evidence has been gathered showing that, for the lowest state of any given symmetry, there is a single nodal hypersurface (up to all permutations) that divides configuration space into only two nodal domains (one positive and one negative) \cite{Ceperley91, Glauser92, Bressanini01, Bressanini05a, Bajdich05, Bressanini05b, Scott07, Mitas06, Mitas08, Bressanini08, Bressanini12}.
Except in some particular cases, electronic or more generally fermionic nodes are poorly understood due to their high dimensionality and complex topology \cite{Ceperley91, Bajdich05}.
The number of systems for which the exact nodes are known analytically is very limited \cite{Klein76, Bressanini05a, Bajdich05, Nodes15}.

The quality of fermion nodes is of prime importance in QMC calculations due to the fermion sign problem, which continues to preclude the application of in principle exact QMC methods to large systems.
The dependence of the DMC energy on the quality of $\PsiT$ is often significant in practice, and is due to the fixed-node approximation which segregates the walkers in regions defined by $\PsiT$ \cite{Ceperley91}.
The fixed-node error is only proportional to the square of the nodal displacement error, but it is uncontrolled and its accuracy difficult to assess \cite{Kwon98, Luchow07a, Luchow07b}.

The DMC method then finds the best energy \emph{for that chosen nodal surface}, providing an upper bound for the ground-state energy. 
The exact ground-state energy is reached only if the nodal surface is exact.
Therefore, one of greatest challenge of QMC methods is to design a well-defined protocol to control the fixed-node error or, equivalently, to be able to build chemical meaningful nodal surfaces for any chemical system.