\documentclass[aip,jcp,preprint,showkeys,column]{revtex4-1}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[french]{babel}
\usepackage{comment}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,hyperref,multirow,amsmath,amssymb,amsfonts,physics,mathpazo,algorithmicx,algorithm,algpseudocode}
\usepackage[normalem]{ulem}
\usepackage[version=4]{mhchem}
\usepackage[compat=1.1.0]{tikz-feynman}

\algnewcommand\algorithmicassert{\texttt{assert}}
\algnewcommand\Assert[1]{\State \algorithmicassert(#1)}
\algrenewcommand{\algorithmiccomment}[1]{$\triangleright$ #1}

\newcommand{\alert}[1]{\textcolor{red}{#1}}
\newcommand{\mika}{Mickaël Véril}
\newcommand{\titou}{Pierre-Fran\c{c}ois Loos}
\newcommand{\lcpq}{Laboratoire de Chimie et Physique Quantiques, Universit\'e de Toulouse, CNRS, UPS, France}

\newcommand{\numsys}[1]{\pmb{\mathbb{{#1}}}}
\newcommand{\evGW}{evGW}	
\newcommand{\qsGW}{qsGW}	
\newcommand{\GOWO}{G$_0$W$_0$}	
\newcommand{\GOW}{G$_0$W}	
\newcommand{\GWO}{GW$_0$}	
\newcommand{\GW}{GW}

\newcommand{\quack}{\texttt{quack}}
\newcommand{\bash}{\texttt{bash}}
\newcommand{\GoDuck}{\texttt{GoDuck}}
\newcommand{\python}{\texttt{python}}
\newcommand{\PyDuck}{\texttt{PyDuck}}
\newcommand{\fortran}{\texttt{fortran90}}
\newcommand{\IntPak}{\texttt{IntPak}}
\newcommand{\MCQC}{\texttt{MCQC}}
\newcommand{\Orca}{\texttt{Orca}}

\newcommand{\Hc}{\hat{H}^\text{c}}
\newcommand{\Fock}{\hat{F}}
\newcommand{\sigc}[1]{\Sigma^\text{c}_{#1}}
\newcommand{\osigc}{\Sigma^\text{c}}
\newcommand{\GOWOSOSEX}{{\GOWO}+SOSEX}		
\newcommand{\GWSOSEX}{{\GW}+SOSEX}	
\newcommand{\GnWn}[1]{G$_{#1}$W$_{#1}$}		
\newcommand{\GOF}{G$_0$F2}
\newcommand{\GF}{GF2}
\newcommand{\nsol}{N_\text{sol}}
\newcommand{\tabref}[1]{(tableau \ref{#1})}
\newcommand{\figref}[1]{(voir figure \ref{#1})}
\renewcommand{\vb}[1]{\bm{#1}}

\newcommand{\x}{\omega-\epsilon_l\pm\Omega}
\newcommand{\y}{\abs{(kl|\rho_m)}}
\newcommand{\re}{\frac{\x}{\y}}

\begin{document}

\title{Mise en évidence de discontinuités dans les surfaces d'énergie potentielle induites par les méthodes \GW}

\author{\mika}
\affiliation{\lcpq}
\email[Corresponding author: ]{mika.veril@wanadoo.fr}

\begin{abstract}
  Cette étude porte sur les méthodes basées sur la fonction de Green, et en particulier l’approximation \GW. 
  Nous démontrons l’existence de discontinuités (non physique) dans les surfaces d'énergie potentielle générées par ces méthodes. 
  Il est donc essentiel de comprendre leurs origines, pour ensuite tenter d'apporter un correctif afin de les éliminer.
  Ces discontinuités sont mise en évidence sur des molécules diatomiques simples, comme \ce{H2}, \ce{HeH+} ou encore \ce{LiF} au niveau auto-cohérent \textit{``eigenvalue-only''} (\evGW).
  L’ensemble des calculs ont été effectués à l'aide de programmes localement développés. 
  Nous avons pu constater que ces discontinuités étaient souvent plus marqué en théorie de la fonctionnelle de la densité (DFT), due a la forte sous-estimation de la différence d’énergie entre la HOMO et la LUMO par ce type de méthodes. 
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Les m\'ethodes bas\'ees sur la fonction de Green, et plus particulièrement l'approximation GW, sont capable de transformer un probl\`eme multi-électronique littéralement impossible \`a r\'esoudre, en un jeu d'\'equations mono-électroniques non linéaire. 
GW prends en compte la corrélation \'electronique gr\^ace a une s\'equence d'\'etapes auto-coh\'erentes, connues sous le nom d'\'equations d'Hedin (voir figure \ref{fig:pentagon} et ci-dessous pour plus amples d\'etails).
D'importantes propri\'et\'es exp\'erimentales peuvent \^etre obtenues \`a partir de la fonction de Green, comme le spectre de photo-\'emission directe ou inverse, ainsi que l'\'energie d'ionisation, l'affinit\'e \'electronique, et bien d'autres observables.
Même si la méthode GW est bien connue pour son application sur les solides, elle est aussi, depuis une dizaine d'années, beaucoup appliquée aux molécules.

Cette étude porte sur les méthodes GW et l’existence de discontinuités dans les surfaces d'énergie potentielle générées par ces méthodes. 
En effet, dans la réalité les surfaces d'énergie potentielle ne comportent bien s\^ur aucune discontinuité. 
Il est donc essentiel de comprendre leurs origines pour ensuite tenter d'apporter un correctif afin de les éliminer.
Plusieurs séries de calculs ont été réalisés.
Une première série de calculs a consisté à faire un ``scan'' de la distance internucl\'eaire ($R$) sur des molécules diatomiques pour trouver des discontinuités dans la surface d’énergie potentielle.
Une seconde série de calculs consistait à faire un calcul GW à la géométrie d'équilibre sur une autre série de molécules (pas nécessairement diatomique), afin d'obtenir la \textit{``self-energy''} $\Sigma$ dans le but de savoir si les géométries d'équilibres sont proches de ces discontinuités.
Une réflexion a aussi été mené pour tenter de trouver des solutions au problème de discontinuités.

Les calculs ont été effectué à l'aide du script {\PyDuck}, que j'ai développé en {\python} autour du logiciel {\quack} de {\titou} du {\lcpq}. 
Ce dernier est constitué d'un script principal {\GoDuck} écrit en {\bash}, du programme de calculs d'intégrales {\IntPak} écrit en {\fortran} et du programme de chimie quantique {\MCQC} également écrit en {\fortran}. 
{\PyDuck} permet d'effectuer des scans, de l'extraction de données et des tracés de graphiques.

%%% Figure 1 %%%
\begin{figure}        
  \includegraphics[width=0.6\linewidth]{fig/fig1}
	\caption{
	  \label{fig:pentagon}
	Pentagone d'Hedin. \cite{Hedin_1970}
	Le chemin rouge correspond au processus GW auto-cohérent qui \'evite le calcul (tr\`es complexe) de la fonction vertex $\Gamma$.}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Théorie}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{L'approximation Hartree-Fock}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
La méthode Hartree-Fock (HF) est une méthode dite mono-\'electronique permettant de résoudre de façon approchée l'équation de Schrödinger indépendante du temps.
Dans l'approximation HF, la fonction d'onde $\Psi_\text{HF}$ s'écrit comme une combinaison linéaire antisymétrique de fonction d'onde monoélectronique appelée spinorbitales $\phi_i(\vb{X})$.
Ceci peut s’écrire élégamment sous la forme d'un déterminant de Slater 
\begin{equation}
  \Psi_\text{HF}=\frac{1}{\sqrt{n!}}
	\begin{vmatrix}
	  \phi_1(\vb{X}_1) & \phi_2(\vb{X}_1) & \cdots  & \phi_n(\vb{X}_1)
		\\
		\phi_1(\vb{X}_2) & \phi_2(\vb{X}_2) & \cdots  & \phi_n(\vb{X}_2)
		\\
		\vdots & \vdots & \ddots &\vdots
		\\
		\phi_1(\vb{X}_n) & \phi_2(\vb{X}_n) & \cdots & \phi_n(\vb{X}_n)
	\end{vmatrix},
\end{equation}
o\`u $n$ est le nombre d'\'electrons du syst\`eme.
Une spinorbitale $\phi_i(\vb{X})=\chi_i(\vb{r})S_i(s)$ est le produit d'une fonction spatiale $\chi_i(\vb{r})$ et d'une fonction de spin
\begin{equation}
  S_i(s)= \begin{cases}
    \alpha, & \text{électron spin-up},
		\\
		\beta, & \text{électron spin-down},
  \end{cases}
\end{equation} 
$\vb{X}$ étant la combinaison du vecteur spatial $\vb{r}$ et du spin $s$.
Les spinorbitales forment une base orthonormale, c'est-\`a-dire $\braket{\phi_i}{\phi_j}=\delta_{ij}$, où 
\begin{equation}
  \delta_{ij}=
   \begin{cases}
     1, & i=j,
      \\
      0,& i \ne j
   \end{cases}
\end{equation}
est le symbole de Kronecker.
On obtient l'opérateur de Fock \cite{SzaboBook} 
\begin{equation}
  \label{eq:Fock}
	\Fock = \Hc + \sum_{j=1}^n \qty( \hat{J}_j-\hat{K}_j ),
\end{equation}
$\Hc$ étant l’Hamiltonien de cœur, c'est à dire la somme de l'op\'erateur énergie cinétique de l’électron et le potentiel d’attraction électron-noyau
\begin{equation}
  \Hc = -\frac{1}{2}\nabla^2 + \sum_{A=1}^{N}\frac{Z_A}{\abs{\vb{r}-\vb{R}_A}},
\end{equation}
o\`u $\hat{J}$ est l'opérateur de Coulomb, $\hat{K}$ l'opérateur d'échange, $N$ le nombre de noyaux, et $\vb{R}_A$ le vecteur position du noyau $A$.
En appliquant l'opérateur de Fock \eqref{eq:Fock} on obtient l'énergie Hartree-Fock $\epsilon_l^\text{HF}$ de l'orbitale $l$, c'est-\`a-dire 
$\epsilon_l^\text{HF} = \ev{{\Fock}}{\phi_l} \equiv \ev{{\Fock}}{l} $.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Fonction de Green}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Généralité}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Comme mentionn\'e précédemment, grâce à un processus auto-cohérent, les m\'ethodes bas\'ees sur la fonction de Green permettent la prise en compte d'une partie de la corrélation électronique. 
La corr\'elation électronique est quantifi\'ee par l'\'energie de correlation, celle-ci étant d\'efinie comme la différence entre l’énergie exacte et l'énergie Hartree-Fock du syst\`eme consid\'er\'e.
Les méthodes de la fonction de Green peuvent se représenter sch\'ematiquement comme sur la figure \ref{fig:pentagon}, \cite{Loos_2018}
o\`u $\Sigma$ est la self-energy, $G$ la fonction de Green, $\Gamma$ la fonction vertex irréductible, $P$ la polarisabilité, $\delta$ la fonction delta de Dirac, $W$ l’int\'eraction écrantée de Coulomb, et enfin $v$ l'int\'eraction de Coulomb (non écrantée).
De plus, il est souvent commode de décomposer la self-energy $\Sigma(\omega) = \Sigma^x\qty(\omega) + \osigc(\omega)$, qui est explicitement d\'ependante de la fr\'equence $\omega$, en une self-energy d'échange statique $\Sigma^x$ et une self-energy de corrélation dynamique $\osigc(\omega)$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{L’approximation GW}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dans l'approximation {\GW} nous utilisons une approximation de la fonction vertex $\Gamma$ (chemin rouge sur la figure \ref{fig:pentagon}):
\begin{equation}
  \label{eq:GammaGW}
	\Gamma(123) = \delta(12) \delta(13).
\end{equation}
o\`u ``1'' (par exemple) est un index composite regroupant espace, spin et temp. 
Ainsi, nous simplifions drastiquement le calcul pour le rendre plus facile à réaliser (ou tout simplement faisable!).
Trois méthodes principales existent: {\GOWO}, {\evGW} et {\qsGW}. 
Nous reportons le lecteur \`a la table \ref{tab:methList} pour de plus amples informations sur ces m\'ethodes.
Ici, nous nous concentrerons sur {\GOWO} et {\evGW}, qui sont de loin les m\'ethodes GW les plus uilis\'ees dans la communaut\'e de la chimie th\'eorique mol\'eculaire.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table*}
  \caption{\label{tab:methList} Liste des méthodes {\GW} les plus populaires.}
  \begin{ruledtabular}
    \begin{tabular}{ccc}
      Acronyme	&	Nom	&	Description	\\
	\hline
	\GOWO	&	GW perturbatif (ou \textit{``one shot''} GW)	&	un seul cycle {\GW} et lin\'earisation			\\
	\evGW	&	\textit{eigenvalue-only} GW 			&	auto-cohérence sur les énergies	orbitalaires		\\
	\qsGW	&	\textit{quasiparticle self-consistent} GW 	&	auto-cohérence sur les orbitales et leur énergies	\\
    \end{tabular}
  \end{ruledtabular}
\end{table*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

L'équation quasiparticule (QP) permet d'obtenir des énergies orbitalaires dites ``quasiparticule'', $\epsilon_l^\text{QP}$, corrig\'ees gr\^ace \`a l'effet de la corr\'elation \'electronique via $\osigc$.
Dans le cas (perturbatif) {\GOWO}, cette \'equation est lin\'earis\'ee, et nous avons
\begin{equation}
  \label{eq:QPO}
	\epsilon_l^\text{QP} = \epsilon_l^\text{HF} + Z_l (\epsilon_l^\text{HF}) \sigc{l}(\epsilon_l^\text{HF}).
\end{equation}
o\`u  le facteur de renormalisation est donn\'e par
\begin{equation}
  \label{eq:Z}
	Z_l(\omega) = \qty[1-\pderivative{\omega}\sigc{l}(\omega)]^{-1},
\end{equation}
avec $0 \le Z_l(\omega) \le 1$.
Quant \`a elle, $\sigc{l}(\omega)$ est obtenue par l'intégrale analytique sur les éléments diagonaux de $\osigc$: \cite{Van_S_2015}
\begin{equation}
  \label{eq:SigR}
	\sigc{l}(\omega)
	= \sum_{m} \qty[ 
    \sum_{i}^\text{occ}
    \frac{\abs{(il|\rho_m)}^2}{\omega - \epsilon_i + \Omega_m}
    + \sum_{a}^\text{virt} 
    \frac{\abs{(al|\rho_m)}^2}{\omega - \epsilon_a - \Omega_m}
    ],
\end{equation}
o\`u ``occ'' et ``virt'' correspondent \`a la somme sur les orbitales occup\'ees et virtuelles, respectivement, et
o\`u $\Omega_m$ est l'énergie d'excitation et $\rho_m$ la densité de transition obtenus par un calcul pr\'eliminaire de type ``random phase approximation'' (RPA).

Dans le cas auto-cohérent {\evGW}, il n'est plus nécessaire d'utiliser le facteur de renormalisation $Z_l$ donn\'e par l'\'equation \eqref{eq:Z} car l'auto-cohérence remplace la lin\'earisation de l'\'equation quasiparticule.
L'energie orbitalaire evGW, $\epsilon_l^\text{QP}$, est obtenu en resolvant l'equation non-lin\'eaire quasiparticule :
\begin{equation}
  \label{eq:QP}
	\omega = \epsilon_l^\text{HF} + \sigc{l}(\omega).
\end{equation} 
De fa\c{c}on pratique, nous iterons de la mani\`ere suivante
\begin{equation}
  \label{eq:QPSC}
	\epsilon_l^\text{QP} = \epsilon_l^\text{HF} + \sigc{l}(\epsilon_l^{\text{QP}-1}),
\end{equation} 
où $\epsilon_l^{\text{QP}-1}$ est l'energie orbitalaire evGW obtenue \`a l'itération précédente.
L'équation \eqref{eq:QP} a bien s\^ur une multitude de solutions. 
Pour savoir quelle est la solution à garder (dite solution quasiparticule), nous suivons g\'en\'eralement la solution ayant le plus grand facteur de renormalisation.
En d'autres termes, nous cherchons parmi les solutions de l’équation \eqref{eq:QP} la valeur de $Z_l$ la plus grande.
Toutes les autres solutions sont ce que l'on appelle des ``satellites''.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Problématique
\label{PB}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L'équation \eqref{eq:SigR} est une combinaison de deux fonctions de la forme $y/x$ avec $y=\abs{(pl|\rho_m)}^2$ et $x=\omega - \epsilon_p \pm \Omega_m$.
$\osigc(\omega)$ n'est donc pas définie pour $\omega - \epsilon_p\pm \Omega_m=0$.
Or dans les \'equations \eqref{eq:QPO} et \eqref{eq:QPSC}, nous évaluons \eqref{eq:SigR} en $\epsilon_l$.
De fa\c{c}on accidentelle, il est tout a fait possible qu'en évaluant la self-energy que le terme $\epsilon_p\pm \Omega_m$ soit tr\`es proche de la valeur $\omega = \epsilon_l$. 
Ceci implique que le dénominateur tend vers $0$ et donc que la fraction tend vers $\pm \infty$. 
Ce pôle de $\osigc$ entraîne l'\'emergence de deux solutions proche en \'energie et donc de discontinuités dans la surface d’énergie potentielle \figref{fig:H2evGW}. 
Ce pose un r\'eel problème car ces discontinuités sont, bien s\^ur, non physique.
Le but de cette étude est donc d'étudier ces discontinuités et les pôles de $\osigc$ qui les engendrent pour comprendre l'envergure du problème et tenter de trouver une approche peu on\'ereuse pour le résoudre.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{PyDuck}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\PyDuck} est un script utilisant {\quack}, un logiciel ``maison'' réalisant des calculs de chimie quantique.
Il se compose d'un script {\python} et d'un modèle de fichier de configuration facilement copiable dans le dossier courant par le biais d'un argument de ligne de commande de {\PyDuck}. 
Ce fichier de configuration est un fichier \texttt{json} (JavaScript Object Notation) \cite{JSONWeb}, format couramment utiliser pour sérialiser des donnés.
Il contient l'ensemble des méthodes prisent en charge par {\quack} avec un commutateur pour activer ou non l'extraction de ces valeurs ainsi que des options additionnelles pour savoir sur quel graphique l'afficher ou quelles orbitales prendre en compte dans le cas d'une propriété orbitalaire.
D'autre options sont présentes concerant l'output ou le mode. 
Il existe actuellement deux modes: i) ``SinglePoint'' où nous pouvons faire des tracés en fonction de $\omega$, et ii) ``Scan'' o\`u les tracés se font en fonction de la distance internucl\'eaire $R$ pour les diatomiques ou de la charge nucl\'eaire $Z$ dans le cas des syst\`emes atomiques. 
Les molécules à plus de deux atomes ne sont pour le moment pas supportées dans ce mode.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
  \includegraphics[width=0.49\linewidth]{fig/eqpH2G0W0}
	\includegraphics[width=0.49\linewidth]{fig/SigCH2G0W0}
	\includegraphics[width=0.49\linewidth]{fig/ZH2G0W0}
	\includegraphics[width=0.49\linewidth]{fig/SigCH2G0W0lumo1}
	 \caption{
	   \label{fig:H2G0W0}
	Calculs {\GOWO}@HF/6-31G sur \ce{H2}.
	Les trois premiers graphiques repr\'esentent $\epsilon_{l}^\text{QP}$, $\sigc{l}(\epsilon_{l}^\text{HF})$ et $Z_l(\epsilon_{l}^\text{HF})$ en fonction de la distance internucl\'eaire $R$.
	Le dernier graphique en bas \`a droite est réalisé vers la discontinuté de la LUMO+1.
	Les solutions de l'\'equation quasiparticule correspondent \`a l'intersection des courbes rouge et noire.
	Nous voyons bien ici qu'il y a 2 solutions tr\`es proche en \'energie.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
  \includegraphics[width=0.49\linewidth]{fig/eqpH2evGW}
	\includegraphics[width=0.49\linewidth]{fig/SigCH2evGW}
	\includegraphics[width=0.49\linewidth]{fig/ZH2evGW}
	\includegraphics[width=0.49\linewidth]{fig/SigCH2evGWlumo1}
	\caption{
	  \label{fig:H2evGW}
	Calculs evGW@HF/6-31G sur \ce{H2}. 
	Le trois premiers graphiques repr\'esentent $\epsilon_{l}^\text{QP}$, $\sigc{l}(\epsilon_{l}^\text{QP})$ et $Z_l(\epsilon_{l}^\text{QP})$ en fonction de la distance internucl\'eaire $R$.
	Le dernier graphique en bas \`a droite est réalisé vers la discontinuté de la LUMO+1.
	Les solutions de l'\'equation quasiparticule correspondent \`a l'intersection des courbes rouge et noire. 
	Nous voyons bien ici qu'il y a 2 solutions tr\`es proche en \'energie.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
  \includegraphics[width=0.49\linewidth]{fig/eqpH2evGWzoom1}
	\includegraphics[width=0.49\linewidth]{fig/eqpH2evGWzoom2}
	 \caption{
	   \label{fig:H2evGWZoom}
	Zoom correspondant aux calculs evGW@HF/6-31G sur \ce{H2}.
	La discontinuit\'e vers $R = 1$ bohr est visible pour la LUMO+2 (graphique de gauche), ainsi que vers $R = 2,25$ bohr pour la LUMO+1 (graphique de droite).}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Résultats}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Scan}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
En faisant varier la distance internucl\'eaire $R$ entre les deux atomes consistuant une mol\'ecule diatomique, nous pouvons identifier, pour différentes molécules, où les discontinuités sont susceptible d'apparaître.
Dans cette section, nous utilisons délibérément une base relativement petite (6-31G) pour pouvoir \'etudier l'ensemble des orbitales.
Dans la section suivante, nous utiliserons une base plus grande.
Nous avons remarqu\'e que l'emplacement des discontinuit\'es est relativement insensible \`a la base.

\paragraph{Cas de \ce{H2}.}
Nous avons trouv\'e plusieurs pôles pour la LUMO+1 et la LUMO+2 (voir figures \ref{fig:H2G0W0} et \ref{fig:H2evGW}).
Dans le premier cas (LUMO+1), nous avons observ\'e deux pôles de $\osigc$ tr\`es proche l'un de l'autre vers $R = 1$ bohr.
Le deuxième cas est plus classique avec un seul pôle au alentour de $R = 2,25$ bohr (voir graphique en haut \`a droite des figures \ref{fig:H2G0W0} et \ref{fig:H2evGW}).
Au niveau evGW, ces pôles entraînent des discontinuités dans $\osigc$ aux mêmes abscisses (voir graphique en haut \`a gauche de la figure \ref{fig:H2evGW}), ce qui entraîne aussi des valeurs du facteur de renormalisation, donn\'e par l'\'equation \eqref{eq:Z}, très petites (voir graphique en bas \`a gauche de la figure \ref{fig:H2evGW}).
Au alentour d'une discontinuit\'e (voir les zooms sur la figure \ref{fig:H2evGWZoom}), le processus auto-coh\'erent peut tomber al\'eatoirement sur l'une ou l'autre des deux solutions de l'\'equation quasiparticule $\omega-\epsilon_l^\text{HF} = \sigc{l}(\omega)$ correspondant \`a l'intersection des courbes rouge et noire sur les graphiques en bas \`a droite dans les figures \ref{fig:H2G0W0} et \ref{fig:H2evGW} (voir section suivante pour plus amples d\'etails).
Ces solutions ayant des \'energies diff\'erentes, ceci implique un saut dans la surface d'\'energie potentielle, donc une discontinuit\'e.
On remarque aussi que plus l'orbitale est haute en énergie plus la discontinuité se présente tôt (c'est-\`a-dire \`a petite distance interatomique $R$). 
Dans le cas de {\GOWO}, nous remarquons que ces pôles n'induisent pas de discontinuités mais une petite bosse appara\^it lorsque $\osigc$ diverge (voir graphique en haut \`a gauche de la figure \ref{fig:H2G0W0}). 
Ceci est d\^u \`a un problème de renormalisation qui n'est pas parfaitement bien réalisé.

\paragraph{Cas de \ce{HeH+} de \ce{LiF}.} 
Pour ces deux mol\'ecules, nous n'avons pas report\'e les graphiques dans ce rapport, mais ceux-ci sont extr\^emement similaire \`a ceux de \ce{H2}.
Pour \ce{HeH+}, nous avons observ\'e deux discontinuités une à $R = 3$ bohr dans la LUMO+2 et une à $R = 4$ bohr dans la LUMO+1.
Dans le case de \ce{LiF}, nous avons uniquement trouv\'e une discontinuité dans la LUMO+1 à $2,5$ bohr.
Ceci montre que \ce{H2} n'est pas un cas isol\'e mais que ce type de probl\`emes appara\^it fréquemment dans les syst\`emes mol\'eculaires.
Cependant, remarquons que, ici, cela ne se produit pas pour les orbitales fronti\`eres (HOMO et LUMO).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{squeezetable}
\begin{table*}
  \caption{\label{tab:gap} Énergie d'ionisation ($\text{IP} = -\epsilon_\text{HOMO}$), affinit\'e \'electronique ($\text{EA} = -\epsilon_\text{LUMO}$), et gap HOMO-LUMO (en eV) obtenus pour certaines molécules ``probl\`ematiques'' étudiées par van Setten et collaborateurs. \cite{Van_S_2015}
	Nos calculs {\GOWO} ont été réalisé avec la base cc-pVDZ en partant d'\'energies orbitalaires obtenues par la méthode Hartree-Fock (HF).
	Pour comparaison, les résultats de van Setten et collaborateurs, réalisés au niveau {\GOWO}@PBE/def2-QZVP, ont aussi été reporté, ainsi que les r\'esultats exp\'erimentaux.}
\begin{ruledtabular}
\begin{tabular}{cccccccccccccccc}
	Mol.		&	\multicolumn{3}{c}{PBE/cc-pVDZ}		&	\multicolumn{3}{c}{HF/cc-pVDZ}	
			&	\multicolumn{3}{c}{{\GOWO}@HF/cc-pVDZ}	&	\multicolumn{3}{c}{{\GOWO}@PBE/def2-QZVP} 
			&	\multicolumn{3}{c}{Exp.}
	\\
	\cline{2-4}	\cline{5-7}	\cline{8-10}	\cline{11-13}	\cline{14-16}
			&	gap	&	EA	&	IP 
			&	gap	&	EA	&	IP 
			&	gap	&	EA	&	IP 
			&	gap	&	EA	&	IP 
			&	gap	&	EA	&	IP 
	\\
	\hline
	\ce{O3}		&	1,71	&	5,59	&	7,30	
			&	12,02	&	1,10	&	13,11	
			&	11,90	&	1,20	&	13,09
			&	9,09	&	2,30	&	11,39
			&	10,63	&	2,10	&	12,73
	\\
	\ce{LiH}	&	2,81	&	1,59	&	4,40 
			&	6,58	&	0,44	&	7,02 
			&	7,16	&	0,31	&	7,47
			&	6,48	&	0,07	&	6,55
			&	7,56	&	0,34	&	7,90
	\\
	\ce{BN}		&	0,19	&	7,11	&	7,30 
			&	8,71	&	2,67	&	11,38 
			&	8,01	&	3,37	&	11,39
			&	7,06	&	3,95	&	11,01 
			&	---	&	3,16	&	---
	\\
	\ce{BeO}	&	1,35	&	4,76	&	6,10 
			&	8,95	&	1,52	&	10,48 
			&	7,43	&	1,94	&	9,37
			&	6,13	&	2,49	&	8,62 
			&	---	&	---	&	10,10	
\end{tabular}
\end{ruledtabular}
\end{table*}
\end{squeezetable}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Single point}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nous avons \'etudi\'e le comportement de $\sigc{\text{LUMO+1}}(\omega)$ et $Z_\text{LUMO+1}(\omega)$ en fonction de la fr\'equence $\omega$ pour les molécules précédemment étudié (\ce{H2}, \ce{HeH+} et \ce{LiF}) au alentour de la distance où se produit une discontinuité. 
Dans cette r\'egion, nous voyons que $Z_\text{LUMO+1}(\omega)$ diminue et $\sigc{\text{LUMO+1}}(\omega)$ présente un pôle \figref{fig:H2evGW}.
La solution de l'équation quasiparticule $\omega-\epsilon_{\text{LUMO+1}}^\text{HF} = \sigc{\text{LUMO+1}}(\omega)$ (équivalent à l'\'equation \eqref{eq:QP}) est symbolisé par l'intersection des courbe de $\osigc$ (courbe noire) et de $\omega-\epsilon^\text{HF}$ (courbe rouge) qui indique les solutions de l'équation quasiparticule.
Ainsi, lorsqu'il y a un pôle, nous avons deux solutions de l'équation quasiparticule \ref{eq:QP} (représenté par les deux intersection du graphique en bas \`a droite de la figure \ref{fig:H2evGW}) avec de petites valeurs de $Z$.
En temps normal, nous avons une solution avec un $Z$ grand ($Z>0,8$) et c'est cette solution qui est retenue; les autres solutions sont appel\'ees satellites.
Dans le cas d'un pôle, deux choix s'offre \`a nous: i) retenir la solution avec le plus grand $Z$ ou ii) retenir la solution la plus proche du niveau de Fermi (qui est souvent d\'efini dans les mol\'ecules comme le point \'equidistant entre la HOMO et la LUMO).
Lorsque $R$ varie, il y a donc à un moment donn\'e un brutal changement de solution d'où la discontinuité quand nous basculons d'un solution à l'autre (voir zooms sur la figure \ref{fig:H2evGWZoom}).

Pour la suite de notre étude, nous nous sommes basés sur les cas problèmatiques rencontrés par van Setten et collaborateurs. \cite{Van_S_2015}
En effet, dans cet article, ils font remarquer que pour certaines molécules, comme \ce{O3}, \ce{LiH}, \ce{BN} et \ce{BeO} \tabref{tab:gap}, au niveau {\GOWO}, qu'ils observent deux solutions de l'\'equation quasiparticule physiquement acceptable pour la HOMO, et ceci \`a la géométrie exp\'erimentale. 
Ceci pose problème car, dans ce cas, nous ne savons pas vraiment quel est la bonne solution de l'équation quasiparticule pour une orbitale qui est directement reli\'ee \`a une grandeur exp\'erimentale (le potentiel d'ionisation).

Concernant les r\'esultats pr\'esents dans le tableau \ref{tab:gap}, une première chose que l'on remarque est que les résultats Kohn-Sham (KS) obtenus avec la fonctionnelle PBE donne, compar\'e a l'exp\'erience, un gap HOMO-LUMO tr\`es petit.
En ce qui concerne nos résultats Hartree-Fock, nous voyons que le gap HOMO-LUMO est bien supérieur au gap HOMO-LUMO PBE, et beaucoup plus proche du gap exp\'erimental.
De plus, les r\'esultats de van Setten au niveau {\GOWO}@PBE/def2-QZVP et les nôtres au niveau {\GOWO}@HF/cc-pVDZ sont assez proche.
Cela montre que la correction {\GOWO}, bien que d\'ependante du point de d\'epart, est efficace et permet de corriger une grande partie des erreurs de la fonctionelle PBE.
Cependant, nos r\'esultats {\GOWO}@HF/cc-pVDZ sont en meilleur accord avec les valeurs exp\'erimentales car notre point de d\'epart est beaucoup plus proche. 

Nous savons que les méthodes GW sont stable au alentour d'environ $3$ fois la valeur du gap ($E_\text{g}$) dans la r\'egion du niveau de Fermi. 
En effet, les premiers pôles de $\osigc(\omega)$ se trouvent à $\omega=\epsilon_\text{HOMO}-\Omega_0$ et $\omega=\epsilon_\text{LUMO}+\Omega_0$.
Or $\Omega_0$ étant la première excitation (HOMO-LUMO), on a approximativement $\Omega_0 \approx E_g = \epsilon_\text{LUMO} - \epsilon_\text{HOMO}$, d'o\`u un intervalle de ``confiance'' de $3 E_\text{g}$.
Ceci implique que nous avons un grand intervalle de valeurs stables autour du niveau de Fermi dans le cas de Hartree-Fock. 
Cependant, nous voyons que dans le cas de PBE, le gap est très petit, donc si la solution de l'équation quasiparticule se trouve en dehors de la limite des $3 E_\text{g}$, des problèmes de discontinuités peuvent apparaître car la solution quasiparticule peut se trouver proche d'un pôle de $\osigc$.
Nous voyons bien sur les graphiques report\'es sur la figure \ref{fig:O3evGW} (et obtenus au niveau {\evGW}@HF/6-31G) cette zone de stabilit\'e autour du niveau de Fermi.
En effet, le facteur de renormalisation $Z_\text{HOMO}(\omega)$ (graphique de droite) est tr\`es proche de l'unit\'e au alentour de la solution quasiparticule (intersection des courbes rouge et noire sur le graphique de gauche).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
  \includegraphics[width=0.49\linewidth]{fig/SigCO3evGW}
	\includegraphics[width=0.49\linewidth]{fig/ZO3evGW}
	\caption{
	  \label{fig:O3evGW}
	Calculs {\evGW}@HF/6-31G pour la HOMO de la mol\'ecule d'ozone \ce{O3} \`a sa géométrie d'équilibre.
	Le graphique de gauche repr\'esente $\sigc{\text{HOMO}}(\omega)$ en fonction de la fr\'equence $\omega$, tandis que le graphique de droite repr\'esente $Z_\text{HOMO}(\omega)$ en fonction de $\omega$.
	Sur le graphique de gauche, la solution de l’équation quasiparticule corresponds a l'intersection des courbes rouge et noire.
	Dans ce cas, cette intersection se situe relativement loin du p\^ole de $\sigc{\text{HOMO}}(\omega)$.}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Les méthodes {\GW} permettent d'obtenir des r\'esultats th\'eoriques en tr\`es bon accord avec les r\'esultats experimentaux.
Cependant, ces discontinuités pr\'esentes dans les surfaces d'énergie potentielle posent un réel problème aux niveaux m\'ethodologique et computationelle.
Nous avons pu montrer que plus l'orbitale est haute plus les discontinuités arrivent tôt (c'est-\`a-dire \`a petit $R$). 
L'autre point important est que la méthode est stable au alentour de 3 fois le gap HOMO-LUMO et comme le gap HOMO-LUMO est sous-estimé au niveau Kohn-Sham et surestimé en Hartree-Fock, ceci implique que Hartree-Fock est une meilleur méthode pour éviter ces discontinuités (tout au moins pour les petits syst\`emes). Pour réellement éliminer ces discontinuités nous pouvons tenter de modifier la fonction $\sigc{l}$ donn\'ee par l'\'equation \eqref{eq:SigR}.
L'id\'ee principale est donc de changer l'expression située dans l'opérateur somme en une expression dont le dénominateur ne diverge pas de fa\c{c}on abrupt.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Remerciements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Je remercie mon encadrant {\titou} pour son aide précieuse tout au long de mon stage, ses explications claires et son soutien, Anthony Scemama pour les bonnes pratiques de fortran qu'il m'a transmis, Arjan Berger pour sa participation active dans mon sujet de recherche. 
Je remercie également tout le personnel du laboratoire pour leur accueil particulièrement Thierry Leininger, Jean-Louis Heully et Eric Colledani. 
Je remercie aussi Léa Brooke et Serigne Sarr mes camarades de promotion et collègues de bureau et aussi mon ancien camarade de M1, actuellement doctorant au laboratoire Mohamed-Amine Bouhammali pour sa gentillesse et son aide ponctuelle.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\bibliographystyle{plainnat}
\bibliography{Rapport}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\appendix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
