%% -*- mode: latex; -*-
\nonstopmode
\documentclass[a4paper]{article}
%% \usepackage[italian]{babel}
\listfiles
%%\usepackage{graphics}
\usepackage{graphicx}
\usepackage{xcolor}
%%\usepackage{amssymb} % e.g. \gtreqqless
\usepackage{amssymb} % e.g. \gtreqqless
\usepackage{amsmath} % e.g. \Re
\usepackage{psfrag}
\usepackage[dvips,all]{xy}
\xyoption{arc}   %% \ellipse; apparently, needs to be loaded explicitly
\xyoption{tile}  % apparently, needs to be loaded explicitly
\usepackage{verbatim}
\usepackage{gblist}
\usepackage{gbshell}
\usepackage[draft]{gbfixme}
\usepackage{ulem}  % overstriking
\usepackage{cancel}  % overstriking in mathmode
% 2015-10-20: apparently, footbib and bibentry can coexist happily
\usepackage[firstcite(page)]{footbib}
\usepackage{bibentry}
\usepackage{units}
\usepackage{adjustbox}
\usepackage{sansmath}

%% \pagestyle{empty}
\parindent=0mm

%% 2017-04-20: to make pdf hyperlinks work: "DVIPS_OPT = -z" in .make.local
%% Moreover, for some obscure reason which I didn't try to clarify,
%% clicking on a page ref link sends me to the wrong page; this seems to
%% get fixed by preceeding every \label with a \phantomsection (defined
%% in hyperref)
\usepackage[debug]{hyperref}
%% \usepackage[ps2pdf]{hyperref}

%% 2017-09-17: NOTE: this MUST come after hyperref, or labels wont be displayed
\usepackage[inline]{showlabels}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\nonstopmode
%%\shorthandoff{"} % UNCOMMENT when both xy and [italian]{babel}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%
\section*{Solid state periodic calculations}
%%%%%%%%%%%%

As already mentioned, in order to build a model system as close as
possible to the one experimentally observed it is important to include
in the calculations the gold surface on which the HAT molecule is
adsorbed.

%%%%%%%%%%%%
\subsection*{Electronic structure calculation for periodic systems}
%%%%%%%%%%%%

Crystalline solids and surfaces can be treated theoretically with the
methods of solid state physics, in which advantage is taken from the
spatial periodicity of the system under study. A solid sample is
considered as a supercell obtained by replication of the
crystallographic primitive cell along the three directions of its
defining edges (figure~\ref{fig:supercell}). By application of the
so--called periodic boundary conditions, the sample becomes in turn
the ``primitive'' cell of an infinitely extended crystal. Thanks to
this definition of the system, it is possible to obtain the
description of the electronic structure of the solid sample free from
any boundary effect and limiting the detailed ab initio atomic
simulation to just the contents of the crystallographic primitive
cell.

\begin{figure}[hb]
\begin{center}
%\fbox%
%                l   b   r   t
{\hspace*{-10mm}\adjustbox{trim=0mm 0mm 0mm 0mm,scale=1.0}{\includegraphics{../../../latex/sss/plot/figures.d/supercell.pdf}}}
\end{center}
\caption{The definition of a solid cristalline sample (the ``large''
  cell marked by bold black edges) as a supercell of the
  crystallographic primitive cell (the yellow cell). The supercell
  edges are integer multiples of the corresponding edges of the
  crystallographic primitive cell:
  {$\vec{A}_{i}=N_{i}\vec{a}_{i}$}. For the supercell in the figure we
  have {$(N_{1},\:N_{2},\:N_{3})=(2,\:4,\:3)$}: of course, for a real
  solid sample, the {$N_{i}$} integers would be of the order of the
  Avogadro number.\label{fig:supercell}}
\end{figure}

In general terms, the procedure to solve the monoelectronic time
independent Schr\"{o}dinger equation for the {$i-$}th electron:

\begin{equation}
  \left(-\frac{\hbar^{2}}{2m_{e}}\nabla^{2}+V\left(\vec{r}\right)\right)\phi\left(\vec{r}\right)=\epsilon\phi_{i}\left(\vec{r}\right)\label{eq:TISE}
\end{equation}

for the case of a crystalline solid where the potential energy
{$V$} is periodic can be summarized as follows.

Since {$V$} is periodic over the Bravais lattice {$\mathcal{R}_{\mu}$}
defining the crystal structure of the solid, it can be expanded in a
Fourier series as:


\begin{equation}
V\left(\vec{r}\right)=\sum_{\vec{G}}V_{\vec{G}}\exp\left(i\vec{G}\cdot\vec{r}\right)\label{eq:V:expansion}
\end{equation}

where {$\vec{G}$} is a vector of the reciprocal lattice
{$\mathcal{R}^{\ast}_{\mu}$} of {$\mathcal{R}_{\mu}$} and:

\begin{equation}
V_{\vec{G}}=\frac{1}{v}\int_{v}V\left(\vec{r}\right)\exp\left(-i\vec{G}\cdot\vec{r}\right)d\vec{r}
\end{equation}

with {$v$} being the volume of the primitive cell of
{$\mathcal{R}_{\mu}$}.

Due to the periodic boundary conditions, the eigenfunction
{$\phi\left(\vec{r}\right)$} is periodic on the Bravais lattice
{$\mathcal{R}_{M}$} defined by the supercell which represents the
solid sample. Note that {$\mathcal{R}_{M}$} is {\it{}commensurate} to
{$\mathcal{R}_{\mu}$}, i.e. the primitive lattice vectors
{$\vec{A}_{i}$} of {$\mathcal{R}_{M}$} are integer multiples of the primitive lattice vectors
{$\vec{a}_{i}$} of {$\mathcal{R}_{\mu}$}:

\begin{equation}
\vec{A}_{i}=N_{i}\vec{a}_{i}~~~~i\in\left\{1,2,3\right\}~~~~N_{i}\in\mathbb{Z}
\end{equation}

Being periodic over {$\mathcal{R}_{M}$}, also {$\phi\left(\vec{r}\right)$}
can be expanded in a Fourier series:

\begin{equation}
\phi\left(\vec{r}\right)=\sum_{\vec{k}}C_{\vec{k}}\exp\left(i\vec{k}\cdot\vec{r}\right)\label{eq:phi:expansion}
\end{equation}

where {$\vec{k}$} is a vector of the reciprocal lattice
{$\mathcal{R}^{\ast}_{M}$} of {$\mathcal{R}_{M}$} and:

\begin{equation}
C_{\vec{k}}=\frac{1}{V}\int_{V}\phi\left(\vec{r}\right)\exp\left(-i\vec{k}\cdot\vec{r}\right)d\vec{r}
\end{equation}

with {$V$} being the volume of the primitive cell of
{$\mathcal{R}_{M}$}, i.e. the volume of the solid sample.

By
plugging~\ref{eq:phi:expansion}~and~\ref{eq:V:expansion}~into~\ref{eq:TISE},
this latter can be converted into an homogeneous system of algebraic
equations for the unknown Fourier coefficients {$C_{\vec{k}}$} and
energy eigenvalues {$\epsilon{}$}:

\begin{equation}
\left(\frac{\hbar}{2m_{e}}k^{2}-\epsilon{}\right)C_{\vec{k}}+\sum_{\vec{G}}V_{\vec{G}}C_{\vec{k}-\vec{G}}=0~~~~\vec{k}\in\mathcal{R}^{\ast}_{M}~~~~\vec{G}\in \mathcal{R}^{\ast}_{\mu}\label{eq:central:equation}
\end{equation}

({$k=\vec{\left|k\right|}$})


Since {$\mathcal{R}_{M}$} is commensurate to {$\mathcal{R}_{\mu}$},
the lattice vectors of {$\mathcal{R}^{\ast}_{\mu}$} are a(n infinite)
subset of those of {$\mathcal{R}^{\ast}_{M}$} and the primitive cell
of {$\mathcal{R}^{\ast}_{\mu}$} contains {$N=N_{1}N_{2}N_{3}$} lattice
points of {$\mathcal{R}^{\ast}_{M}$}.

Even if this would not be mandatory in principle, it is most
convenient to take the first Brillouin zone, {$\mathcal{B_{Z}}$}, as
the the primitive cell of {$\mathcal{R}^{\ast}_{\mu}$} among all the
infinite possible choices.

By taking advantage of the fact that {$\mathcal{R}_{M}$} is
commensurate to {$\mathcal{R}_{\mu}$} it is possible to partition the
system~\ref{eq:central:equation} into {$N$} independent subsystems
indexd by the {$N$} lattice points of {$\mathcal{R}^{\ast}_{M}$}
inside {$\mathcal{B_{Z}}$}:

\begin{equation}
\left(\frac{\hbar}{2m_{e}}\left|\right.\vec{k}-\vec{G}\left.\right|^{2}-\epsilon{}\right)C_{\vec{k}-\vec{G}}+\sum_{\vec{G}^{\prime}}V_{\vec{G}^{\prime}}C_{\vec{k}-\vec{G}^{\prime}}=0~~~~\vec{k}\in\mathcal{B_{Z}}~~~~\vec{G},\:\vec{G}^{\prime}\in\mathcal{R}^{\ast}_{\mu}\label{eq:subsystem:k}
\end{equation}


The original problem has thus been converted into {$N$} algebraic
systems, one for each {$\vec{k}$} point inside {$\mathcal{B_{Z}}$}.

For a fixed {$\vec{k}$} point inside {$\mathcal{B_{Z}}$}, the solution
of the corresponding system will give an infinite (in principle) set
of eigenvalues {$\epsilon{}$} and associated eigenfunctions
{$\phi\left(\vec{r}\right)$} (i.e. a set of Fourier coefficients
{$C$}). The eigenvalues sorted in ascending order can be indexed by a
single integer {$n$} (called {\it{}band index}) and a single solution
of~\ref{eq:TISE} is thus completely indexed by the {$\vec{k}$} point
inside {$\mathcal{B_{Z}}$} and the band index {$n$}:

\begin{equation}
\epsilon{}_{\vec{k}n}~~~~\phi_{\vec{k}n}\left(\vec{r}\right)\nonumber{}
\end{equation}

%%%%%%%%%%%%
\subsection*{Limiting the number of
  \raisebox{1pt}{{$\stackrel{\mathbf{\rightarrow}}{\mathbf{G}}$}} vectors to be
  considered}
%%%%%%%%%%%%

The system~\ref{eq:subsystem:k} has as many equations as {$\vec{G}$}
vectors, i.e. an infinite number. To make it computationally
affordable only the finite number of {$\vec{G}$} vectors inside an
energy isosphere in {$\vec{G}$} space is usually considered:

\begin{equation}
\vec{G}:~~~~\frac{\hbar{}\left|\right.\vec{k}-\vec{G}\left.\right|^{2}}{2m_{e}}\le{}E_{\textrm{cutoff}}\label{eq:cutoff}
\end{equation}

where {$E_{\textrm{cutoff}}$} is a predefined cutoff energy value.

Of course, the precision of the calculation will be related to 
{$E_{\textrm{cutoff}}$}, larger values implying more {$\vec{G}$}
vectors considered and hence a higher precision. On the other hand,
increasing the cutoff value makes the computation more resource
demanding, so a balance has to be found.

%%%%%%%%%%%%
\subsection*{Pseudopotentials}
%%%%%%%%%%%%

Another crucial point regarding the computational implementation of
the solution of system~\ref{eq:subsystem:k} is related to the
potential energy {$V\left(\vec{r}\right)$}.

In the frame of density functional theory, {$V$} is written as:

\begin{equation}
V=V_{\textrm{nuc}}+V_{H}+V_{xc}\label{eq:total:potential}
\end{equation}

where:
\begin{gblist}{XXXXX}
\item[{$V_{\textrm{nuc}}$}] is the nuclear potential acting on the
  electron
\item[{$V_{H}$}]  is the Hartree potential
\item[{$V_{xc}$}] is the exchange--correlation potential
\end{gblist}

Due to the {$-\infty$} divergence of {$V_{\textrm{nuc}}$} at the
nuclear positions ({$V_{\textrm{nuc}}\approx-Ze^{2}/r$}, {$Z$}=atomic
number) the atomic core wavefunctions display wild high frequency
oscillations when approaching the nucleus and these would require an
enormous number of plane waves to be properly represented by the
Fourier expansion.

Roughly, for a function oscillating with a length scale of
{$\delta$}, a Fourier expansion with wave vector components up to
{$2\pi/\delta$} is required.

As an example, the {$1s$} atomic wave function of the {$C$} atom has
an associated length scale of approximately
{$\delta\approx\unit[0.1]{a.u.}$}. For {$C$} in the diamond polymorph
(volume of the primitive cell: {$v=\unit[76.5]{a.u.}\implies$}~volume of
the reciprocal primitive cell ({$=$}~volume per {$\vec{G}$} vector in
reciprocal space): {$8\pi^{3}/v$}) the number of required
{$\vec{G}$} vectors in the sum of~\ref{eq:subsystem:k} would be:

\begin{equation}
\frac{\textrm{volume~of~a~sphere~of~radius}~\dfrac{2\pi}{\delta}~\textrm{in~reciprocal~space}}{\textrm{volume~per~{$\vec{G}$}~vector}}=\frac{\dfrac{4}{3}\pi\left(\dfrac{2\pi}{\delta}\right)^{3}}{\dfrac{8\pi^{3}}{v}}\approx320000\nonumber{}
\end{equation}

\begin{comment}
% octave start
delta=0.1;
v=76.5;
((4/3)*pi*(2*pi/delta)^3)/(8*pi^3/v)
% octave end
ans = 320442.4506661589
\end{comment}

In order to cut the number of plane waves in the Fourier expansion of
the potential to a manageable value, the so--called pseudopotential
technology has been developed.

Basically, the pseudopotential is a smooth effective nuclear potential
able to reproduce the effect of the nucleus plus core electrons on the
valence electrons.

The idea is that most of the properties one is interested in when doing
an electronic structure calculation depend almost exclusively upon the
behavior of the outermost (valence) electrons of the atoms,
i.e. electrons far from the nucleus.

As far as this is assumed to be the case, then replacing real atomic
orbitals with {\it{}pseudo-orbitals} matching the behavior of the real
ones as closely as possible outside a given distance from the nucleus
should not have any sizeable effect on the calculated properties.

On the other hand, if one is able to make these pseudo-orbitals behave
smoothly inside a given core region (while maintaing their
``faithful'' behavior outside it), the number of Fourier terms needed
to represent them can be dramatically reduced, thus making the
system~\ref{eq:subsystem:k} computationally approachable.

As said above, the goal of smoothing the one electron orbitals inside
the core region while keeping their original behavior outside it is
achieved by replacing the true nuclear potential {$V_{\textrm{nuc}}$}
with a properly designed pseudopotential.

The generation of atomic pseudopotentials is far from trivial and this
technology has been (and is still being actively) developed over a
long period of time. To better clarify the concept of a
pseudopotential, figure~\ref{fig:pseudopotential} shows the
all-electron and pseudo counterparts for the nuclear potential and {$4s$}
orbital of a {$Ti$} atom. The diverging behavior of the 
potential felt by the {$Ti$} {$4s$} electron is replaced by a finite
pseudopotential: the consequence is that the oscillations of the
atomic orbital in the core region are eliminated, while maintaining
the original behavior at distances {$\ge\approx\unit[2]{bohr}$} from
the nucleus.

\begin{figure}[h]
\begin{center}
%% \fbox%
%% NOTE: this needs \usepackage{graphicx}
%                l   b   r   t
{\adjustbox{trim=30mm 0mm 30mm 0mm,scale=1.0}{\sf\sansmath\includegraphics{./figures.d/pseudo-potential.pdf}}}
\end{center}
\caption{Lower panel: (local part of) the atomic potential for a
  {$4s$} electron in the {$Ti$} isolated atom; dotted line is the true nuclear
  potential; full line is the pseudopotential. Upper panel: radial
  part of the {$4s$} orbital; dotted line is the undumped behavior;
  full line is  the pseudo orbital behavior.\label{fig:pseudopotential}}
\end{figure}


%%%%%%%%%%%%
\subsection*{Practical implementation of the calculations}
%%%%%%%%%%%%

We have performed the periodic calculations needed for this thesis
using the Quantum Espresso (QE) [{\bf{}FIXME}] suite of codes. This is
an actively maintained and internationally recognized software for the
practical implementation of the theoretical aspects outlined above.

The {$V_{\textrm{nuc}}$} contribution to the total potential {$V$} (eq.~\ref{eq:total:potential}) is
built with carefully designed atomic pseudopotentials, which are given
as input quantities. The remaining two contributions, {$V_{H}$} and
{$V_{xc}$} depend on the electron density {$n\left(\vec{r}\right)$}:

\begin{eqnarray}
V_{H}&=&V_{H}\left(n\left(\vec{r}\right)\right)\\
V_{xc}&=&V_{xc}\left(n\left(\vec{r}\right)\right)\\
\nonumber
\end{eqnarray}

In turn, the electron density can be calculated only after the one electron
orbitals have been obtained by the solution of the set of
monoelectronic Schr\"{o}dinger equations~\ref{eq:TISE}:

\begin{equation}
n\left(\vec{r}\right)=\sum_{i}^{N}\left|\phi_{i}\left(\vec{r}\right)\right|^{2}
\end{equation}

For this reason, as is usually the case in general, the set of
monoelectronic Schr\"{o}dinger equations~\ref{eq:TISE} is 
solved self--consistently for a given set of {$\vec{k}$} points given
as input. For each {$\vec{k}$} point:

\begin{gblist}{XX}
\item[\textbullet] the {$V_{\textrm{nuc}}$} potential is constructed with
  the atomic pseudopotentials
\item[\textbullet] a starting guessed electron density is built up
\item[\textbullet] {$V_{H}\left(n\left(\vec{r}\right)\right)$} and
  {$V_{xc}\left(n\left(\vec{r}\right)\right)$} are calculated
\item[\textbullet] the monoelectronic orbitals
  {$\phi_{i}\left(\vec{r}\right)$} are obtained by converting
  eq.~\ref{eq:TISE} into the corresponding Fourier space
  version~\ref{eq:subsystem:k}. A crucial computational point in this
  regard is the use of Fast Fourier Transform techniques for the
  generation of the {$V_{\vec{G}^{\prime}}$} Fourier coefficients on a
  grid defined by the cutoff energy condition~\ref{eq:cutoff}.
\item[\textbullet] a new electron density can now be calculated:
{$n\left(\vec{r}\right)=\sum_{i}^{N}\left|\phi_{i}\left(\vec{r}\right)\right|^{2}$}
\item[\textbullet] the process is iterated until the electron density
  at step {$j$} differs from that at step {$j-1$} by less than a
  predefined convergence threshold 
\end{gblist}

The whole process is represented in the flowchart shown in
figure~\ref{fig:flowchart}.


\begin{figure}[h]
\ifpdf
\begin{center}
%%\fbox%
%%                  L   B   R   U
{\adjustbox{trim=0mm 0mm 0mm 0mm,scale=1.0}{\sf\includegraphics{figures.d/flowchart.pdf}}}
\end{center}
\else
\begin{center}
\fbox{{\bf{}FIXME}}
\end{center}
\fi
\caption{The self consistent field procedure as implemented in the
  Quantum Espresso code.\label{fig:flowchart}}
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
\endinput
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Local Variables:
eval: (load-file "./.macros.el")
End:
