~b-ly/+junk/masterThesis

 First-principles studies, also known as \textit{ab initio}, do not require the knowledge of experimental data. Studies are performed ``from scratch'' thanks to ``first-principles''. Only fundamental constants are needed to perform the analysis of a system. Namely, \abinitio\ calculations use:

Last, \siesta\ implementation of the theory is presented. It is based on the usage of a local basis set. \siesta\ can be seen as a tight-binding method which takes into account the presence of all the neighbours. 

\siesta\ wavefunctions are described on the basis of atomic orbitals (LCAO).

  \input{pstricksImg/exHole}

  \caption[Electron hole]{Due to its intrinsic properties, electrons create a depletion around themselves.}

    \caption[Pseudopotential of silicon]{Beyond a given cutoff radius (1.9 in this case), all pseudopotentials are equal to the ionic potential which is represented by dashed line. The pseudos were generated using the atom program delivered with \siesta\ following the scheme of \textcite{Troullier1993}. These pseudos are used throughout this work.}

\section{\siesta\ implementation: atomic orbitals}

The \siesta\ project (Spanish Initiative for Electronic Simulations with Thousands of Atoms) was born in 1996 to provide an order-$N$ simulation package within the self-consistent DFT framework \cite{siestaWebsite}. 

The aim of the current section is to explain and detail how \siesta\ implements the DFT framework and to understand the differences between linear and cubic scaling. Focus is put upon the specificities of \siesta: the usage of atomic orbitals as basis set. A complete description of the \siesta\ method can be found in the literature \cite{PhysRevB.53.R10441,siesta2002}, a more recent article describes the recent developments of the project \cite{Artacho2008}. Initial articles explaining how the basis set is generated are also worth reading \cite{PhysRevB.64.235111,Artacho1999,PhysRevB.66.205101}.

The current section is based on those articles and on workshops organised by the \siesta\ project.

Then, \siesta's bases are detailed along with other specific choices made by the project. A comparison with planewaves id also given.

Linear scaling with respect to the number of atoms was one the objective of \siesta\ and the choices made by the project clearly reflect this goal. The key to linear scaling is the explicit use of locality. Locality means that perturbations which are ``far'' away from a region of the system do not significantly influence properties of that region. While this ``near-sightedness'' principle \cite{Kohn1996} seems straightforward for some systems, it is less evident in the case of metallic compounds. 

Solving the electronic problem can be roughly decomposed into two steps. First the Hamiltonian needs to be built. Second, it must be solved. Each step is done with specific algorithms and do not necessarily have the same type of scaling. A large amount of researchers have worked on improving the algorithms to solve the Hamiltonian. On the other hand, \siesta\ focuses on the first step which afterwards permits to achieve the order-$N$ goal. Once \siesta\ has built the Hamiltonian, it is up to the user to choose the way it is solved:

with $c_1 > c_3$. Figure \ref{fig:cpuTime} illustrates the evolution of the computation time (CPU) with respect to the number of atoms in the system. For small systems, it is more interesting to use diagonalization method whereas very large systems clearly benefit from linear scaling.

  \caption[Linear and cubic scaling of \siesta\ algorithms]{The choice for the solving method depends on the size of the system. For small systems, it is advantageous to use ``classical'' diagonalization technique while for larger systems the computation time is strongly decreased when using linear-scaling algorithms.}

In order to build an Hamiltonian suited for order-$N$ algorithms, a large amount of strategies have been proposed. In his PhD thesis, \textcite{Haynes1998} presents Bessel functions in overlapping spheres as candidates. Another approach is the use of systematic basis set consisting of localized functions on the points of a three-dimensional grid. Such functions are known as ``blips'', they are centered on the points of a grid, and vanish exactly outside a limited domain surrounding each grid point. This technique was proposed by \textcite{PhysRevB.55.13485} and has been implemented in the {\sc Conquest} code \cite{conquest}.  \textcite{PhysRevB.62.1713} have chosen a method based on a real space grid with a finite- difference method. Even a plane-wave approach has been proposed and implemented to provide order-$N$ scaling by making well controlled approximations \cite{onetep}. The last approach, which is the one chosen by \siesta\ developers, is to take advantage of atomic orbitals. 

Different techniques can be used in order to represent the orbitals of the basis set. A popular choice is the usage of Gaussian type orbitals (GTO). Another possibility is the usage of Slater type orbitals. The \siesta\ project is implemented with Numerical Atomic Orbitals (NAO). \\ The purpose of the next section is to detail this approach. 

  In order to build the Hamiltonian, \siesta\ uses numerical atomic orbitals (NAO) as a basis set. Molecular orbitals representing the bonding between two atoms are described through a Linear Combination of Atomic Orbitals (LCAO). 

NAO are the numerical solutions of the KS Hamiltonian for the isolated atom with the same approximations (exchange-correlation and pseudopotential) as the ones that are used for the condensed system. When pseudopotentials are used the basis set is composed of pseudo-atomic orbitals (PAO).

 

 

  The main advantage of atomic orbitals is their efficiency, only a small number of orbitals are needed to reach a good precision. Other advantages include direct interpretation of the results as concepts of orbitals are inherent in chemistry. On the other hand, orbital bases present the disadvantage of lacking a systematic method for optimal convergence.

 

  \siesta\ is able to generate a basis set given a small number of parameters but it is not exclusive in that users can perform simulations with their own orbitals and basis set. Two conditions are to be met though

  If the user does not provide its own basis, \siesta\ will generate it. For this purpose, the two most important parameters are:

    \item[The range] of the orbitals is controlled by a single parameter expressing the energy shift (increase) induced by limiting its range while keeping its norm. 

  These are not the only parameters through which the generation of the basis is controlled. For example, it is possible to control the way multiple-$\zeta$ are generated. The \siesta\ manual provides the user all the required information in this regard \cite{siestaMan}.\\

  \subsubsection{Range of the orbitals}

\siesta\ implements an universal method to obtain finite-range orbitals. Beyond a cutoff radius, orbitals goes to zero. This is performed by implementing the method of \textcite{PhysRevB.40.3979}: single-$\zeta$ orbitals are numerical solutions of:

%This procedure is performed in such a way that the norm of the integral is conserved. This results in an energy shift (increase). By defining this shift in energy, the cutoff radius is computed. 

\begin{equation}

  \left(

    -\frac{1}{2r} \ddtot{}{r} r 

    + \frac{l \left(l+1 \right)}{2r^2}

    + V_l(r)

  \right) \phi(r)

  = \left(\e_l + \Delta E_{\text{PAO}} \right)\phi(r),

  \label{eq:deltaPAO}

\end{equation}

in which $V_l(r)$ is the atomic pseudopotential and $\Delta E_{\text{PAO}}$ is the energy shift  such that $\phi(r)=0$ for all $r>R^l_c$.

The larger the allowed shift, the more confined are the orbitals. A small energy shift thus allows the orbitals to be more extended. This method which is universal in that it can be applied to any element without restriction, is illustrated in \figref{fig:cutoffRadiusEnergyShift}. 

While it is systematic or universal, the method does not guarantee to provide the most suited orbitals for specific problem. It is up to the user to test the basis and make sure that the studied system is correctly described. 

 

\begin{figure}

  \centering

  \input{pstricksImg/energyShit}

  \caption[Range of atomic orbitals controlled by the energy shift ($\Delta E_{\text{PAO}}$)]{\siesta\ is able to generate finite-range orbitals by specifying a single parameter: the energy shift ($\Delta E_{\text{PAO}}$). It expresses the energy increase implied by the confinement of the orbitals. The case of silicon is illustrated (only the radial part of the orbital is plotted). {\color{gray} $R_c(500)$} and $R_c(50)$ are to the cutoff radii corresponding to {\color{gray} $\Delta E_{\text{PAO}} = 500$} and $\Delta E_{\text{PAO}} = 50$ respectively.}

  \label{fig:cutoffRadiusEnergyShift}

\end{figure}

 

 

 

  \siesta\ follows the nomenclature of Quantum Chemistry to name the different options regarding the size of basis sets. 

    \item[A Multiple-$\zeta$] basis follows the same rule as single-$\zeta$ regarding which angular momentum channels are described, but extends the number of radial functions describing them. In the case of a double-$\zeta$ basis, two radial functions are used per valence-orbital. As each radial function attached to a valence-orbital has a different spatial expansion, multiple-$\zeta$ offer radial flexibility.

  Each ``add-on'' permits to better describe the electronic density. By increasing the basis size, the precision is improved. As the size of the basis can be increased at will, there exists no intrinsic limitation for LCAO with respect to quality. However, increasing the basis size can not be performed systematically in that there are different ways to increase the basis size, there is no single parameter like is the case with planewaves for which the cutoff energy can be changed at will in a systematic manner.

  However, one must always check whether the quality improvements are worth the computation time. 

   \siesta\ implementation of the split scheme is illustrated in \figref{fig:splitScheme}. The procedure can be divided into three steps. First a function reproducing the tail of the first-$\zeta$ orbital beyond a given radius, $R_{\text{DZ}}$ and, which smoothly goes to the origin is produced:

   in which, the two parameters ($a$ and $b$) permit to ensure that this function and its first derivative are continuous at $R_{\text{DZ}}$. $l$ refers to the ``$l$'' quantum number.

     \bar{R}^{2\zeta} = R^{1\zeta} - r^l~\left(a - b r^2 \right)

   Finally, the latter function ($\bar{R}^{2\zeta}$) is normalized to obtain the double-$\zeta$ orbital. The value of the radius $R_{\text{DZ}}$ is controlled by the splitNorm parameter inside the \siesta\ code. This parameter thus controls the spatial expansion of the multiple-$\zeta$ orbitals which is always smaller than the one of the first-$\zeta$ orbitals. 

     \caption[Generation of multiple-$\zeta$ bases]{Multiple-$\zeta$ bases have radial flexibility and can be generated by the split scheme using numerical atomic-orbitals. Generation of double-$\zeta$ is illustrated in the case of silicon. The cutoff radius, $R_{\text{DZ}}$ is controlled by the PAO.SplitNorm parameter and is always smaller than $R_c$.}

  Polarization can be achieved in different ways. To polarize a shell with momentum $l$, one can simply add a shell $l+1$. However, this technique results in orbitals that are too extended. Therefore, \siesta\ adopts another method. Numerical orbitals, resulting from the actual polarization of the pseudo-atomic orbital in the presence of a small electric field, are solved exactly with DFT. This leads to $l+1$ orbitals by comparison with first order perturbation theory. \\

  If a double-$\zeta$ basis is used, each type of orbital is doubled: two $s$, two $p_x$, two $p_y$ and two $p_z$. Thus, double-$\zeta$ basis set gives radial flexibility by adding one more radial function within the same angular momentum shell. In this case, the number of orbitals is equal to height. 

  \caption[Size of \siesta\ bases in the case of silicon]{The size of the basis set can be determined through the usage of multiple-$\zeta$ (radial flexibility) and polarized (angular flexibility) orbitals. The accuracy of the results is improved as the basis size is increase, but so does the computation time. The case of silicon is illustrated, a polarized single-$\zeta$ basis will thus have $4+5=9$ orbitals per silicon atom, while for a polarized double-$\zeta$ this number grows to $8+5=13$.}

In section \ref{sec:bsSiBulk}, the effect of the basis size on the bandstructure of bulk silicon is studied. The results are summarized on figure \ref{fig:bsComparison} and table \ref{tab:basisComparison}. Polarization and multiple-$\zeta$ extend the basis in different ways and have thus a different influence on the accuracy. In the case of bulk silicon, using a polarized basis is more important than using a double-$\zeta$ one. This is particularly seen when the bandstructure of silicon is studied.

Using planewaves (PW) to construct the basis set constitutes another popular way to implement the DFT, and a lot of codes have chosen this route. To name only one, the {\sc Abinit} package is developed here  at ``Université catholique de Louvain'' and uses planewaves as basis functions. planewaves present the advantage that a systematic convergence method with respect to the basis can be applied easily. The number of PW composing the basis is controlled by an energy cutoff which represents the maximum kinetic energy a plane wave can have and for which no aliasing effects are observed during Fourrier transforms. While the basis size can be large, they allow simple evaluation of matrix elements and easy parallelization. PW methods present the disadvantage of scaling with the volume of the periodic cell and not directly with the number of atoms. Indeed, the number of planewaves for a given cutoff scales with the volume. Therefore, computations for systems which comprise a large amount of free volume may become very expensive. Such systems include low-density porous structures, molecules or 2D surfaces which is our interest here. Indeed, the latter have to be put into a sufficiently large simulation box in order to prevent interactions between repeated/translated replicas of that box. On the other hand, atom-centered bases scale directly with the number of atoms and the size of the basis (see section \ref{subsec:paoBasis} and table \ref{tab:siliconBasisSets}). In this case, adding free volume inside the unit cell barely costs anything as the basis functions have a limited range and will not describe this area. If one could have the guarantee for a good atom-centered basis, that is, a basis which offers good variational freedom in the regions where it is most needed, i.e. where bonds are created between atoms, then computations would be very efficient. The problem is to find that ``good'' basis, and this is precisely not an easy task. \siesta\ offers a systematic way to generate basis functions based on two main parameters: the basis size and the energy shift implied by its finite range. However, remember that the quality of the basis is not guaranteed. This method does not provide the optimized basis for all atoms and systems. Therefore, one should always check the basis.\\

Within \siesta, the matrix elements are evaluated using planewaves on a real-space grid. Hence, \siesta\ has also a concept of energy cutoff (controlled by the MeshCutoff parameter). 

  \caption[Grid cutoff in real space]{The grid cutoff determines the maximum kinetic energy of the planewaves that can be represented in this grid without aliasing.}

The grid cutoff from \siesta\ should not be confused with the cutoff energy of PWs. Indeed, in the case of a PW basis, electronic wavefunctions are described as a linear combination of PWs and those PWs are projected on a real-space grid. On the other hand, \siesta\ uses a real grid only to represent the charge density.

%\subsection{Obtain STM images within \siesta}