\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{\mika}{Mickaël Véril}
\newcommand{\titou}{Pierre-Fran{\c c}ois Loos}

\newcommand{\lcpq}{Laboratoire de Chimie et Physique Quantiques}
%\newcommand{LCPQCNRS}{{LCPQ}, 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{fortran}}
\newcommand{\IntPack}{\texttt{IntPak}}
\newcommand{\MCQC}{\texttt{MCQC}}
\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{\Orca}{\texttt{Orca} \cite{Neese_2012}}
\newcommand{\figref}[1]{(voir figure \ref{#1})}
\renewcommand{\vb}[1]{\bm{#1}}
\begin{document}
\title{Mise en évidence de discontinuités dans les surfaces d'énergie potentielle induites par les méthodes auto-cohérente \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. \cite{Loos_2018} Ces discontinuités sont mise en évidence sur des molécules diatomiques simples, comme H2, \ce{HeH+} ou encore \ce{LiF} au niveau perturbatif GW (G0W0) et auto-cohérent \evGW. \cite{Bruneval_2016} L’ensemble des calculs ont été effectués à l'aide de programmes localement développés. Une réflexion est aussi menée pour tenter de trouver une solution aux problèmes de discontinuités. En particulier, nous proposons une méthode simple afin de les atténuer. Nous avons pu constater que ces discontinuités était 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éthode. 
\end{abstract}
\maketitle

\section{Introduction}
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 aucune discontinuité. Il est donc essentiel de comprendre leurs origines  ensuite tenter d'apporter un correctif afin d'éliminer ces discontinuités. Plusieurs séries de calculs ont été réalisés ;
une première série de calculs a consisté à faire un scan sur des molécules diatomiques pour trouver des discontinuités dans la surface d’énergie potentielle en fonction de la distance.
Une seconde consistait à faire un calcul GW à la géométrie d'équilibre sur une autre série de molécules (pas nécessairement diatomique), pour tracer la self-energy dans le but de savoir si les géométries d'équilibres sont proches des 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 {\IntPack} é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.
\section{Théorie}
\subsection{Hartree-Fock}
La méthode Hartree-Fock (HF) est une méthode 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(\vb{X})$.
Ceci peut s’écrire élégamment sous la forme d'un déterminant de Slater avec pour $n$ le nombre d'électrons.
\begin{align}
\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}
\intertext{Une spinorbitale $\phi_i$ est le produit de la fonction spatiale $\chi_i$ et de la fonction de spin $S_i$, avec $\vb{X}$ étant la combinaison du vecteur spatial $\vb{r}$ et du spin $s$}
\phi_i(\vb{X})=\chi_i(\vb{r})S_i(s)
S_i(s)= \begin{cases}
      \alpha, & \text{électrons spin-up} 
      \\
      \beta, & \text{électrons spin-down}
\end{cases}
\intertext{Les spinorbitales forment une base orthonormale}
\braket{\phi_i}{\phi_j}=\delta_{ij}
\end{align}
Où $\delta_{ij}$ est le symbole de Kronecker
\begin{equation}
\delta_{ij}=
   \begin{cases}
      1, & i=j 
      \\
      0,& i \ne j
   \end{cases}
\end{equation}
On obtient l'opérateur de Fock 
\begin{align}
\label{eq:Fock}
\Fock=\Hc+\sum_{j} \qty( \hat{J}_j-\hat{K}_j )
\\
\intertext{$\Hc$ étant l’hamiltonien de cœur, c'est à dire la somme de l'énergie cinétique de l’électron et du potentiel d’attraction électron-noyaux}
\Hc=-\frac{1}{2}\nabla^2+\sum_{m=1}^{N}\frac{Z_m}{\abs{\vb{r}-\vb{R}_m}}
\end{align}
on à
$\hat{J}$ l'opérateur de Coulomb
$\hat{K}$ l'opérateur d'échange, $N$ le nombre de noyaux, $\vb{r}$ la position du noyaux dans l'espace et $\vb{R}_m$ la position du noyaux $m$ dans l'espace.\\
En appliquant l'opérateur de Fock \eqref{eq:Fock} on obtient l'énergie Hartree-Fock $\epsilon_l^\text{HF}$ de l'orbitale $l$
\begin{equation}
\epsilon_l^\text{HF}=\ev{{\Fock}}{\phi_l}
\end{equation}
\subsection{Fonction de Green}
\subsubsection{Généralité}
Grâce à un processus auto-cohérent, la fonction de Green permet la prise en compte d'une partie de la corrélation électronique. Celle-ci étant la différence entre l’énergie exacte et l'énergie Hartree-Fock.
\\
La méthode peut se représenter graphiquement \figref{fig:pentagon}
%%% figure 1 %%%
\begin{figure}        \includegraphics[width=0.4\linewidth]{fig/fig1}
        \caption{
        \label{fig:pentagon}
        Pentagone d'Hedin. \cite{Hedin_1970}
        Le chemin rouge montre le processus GW auto-cohérent qui simplifie le calcul de la fonction vertex $\Gamma$.}
\end{figure}
%%%     %%%
\begin{subequations}
  \begin{align}
    \label{eq:Sig}
	& \Sigma(12) = i \int G(13) W(14) \Gamma(324) d(34),
    \\
	\label{eq:G}
	& G(12) = G_\text{HF}(12) + \int G_\text{HF}(13) \Sigma(34) G(42) d(34),
	\\
	\label{eq:Gamma}
	& \Gamma(123) = \delta(12) \delta(13) 
	\notag 
	\\
	& \qquad \qquad		+ \int \fdv{\Sigma(12)}{G(45)} G(46) G(75) \Gamma(673) d(4567),
	\\
	\label{eq:P}
	& P(12) = - i \int G(13) \Gamma(324) G(41) d(34),
	\\
	\label{eq:W}
	& W(12) = v(12) + \int v(13) P(34) W(42) d(34),
  \end{align}
\end{subequations} \cite{Loos_2018}
avec $\Sigma$, la self energy ;
$G$ la fonction de Green ;
$\Gamma$ la fonction vertex irréductible ;
$P$ la polarisabilité
$\delta$ étant la fonction delta de Dirac ;
$W$ l’interaction écrantée de Coulomb.
et $v$ est l'interaction de Coulomb;
Chaque chiffre représente un électron avec ses coordonnées spatiale, son spin et sa coordonnée temporelle.
\\
On peut décomposer la self-energy $\Sigma$ en une énergie d'échange $\Sigma^x$ et l'énergie de corrélation $\osigc$
\begin{equation}
\Sigma=\Sigma^x+\osigc
\end{equation}
\subsubsection{L’approximation GW}
Dans l'approximation {\GW} nous utilisons une approximation de $\Gamma$ à la place de l'équation \eqref{eq:Gamma} pour simplifier cette étape de calculs :
\begin{equation}
        \label{eq:GammaGW}
        \Gamma(123) = \delta(12) \delta(13)
\end{equation}
on obtient l'équation \eqref{eq:GammaGW} en ne prenant uniquement que le premier terme de la partie droite de l'équation {\eqref{eq:Gamma}}.
Ainsi nous simplifions le calcul pour le rendre plus facile à réaliser \figref{fig:pentagon}.
Même si cette méthode est bien connue pour son application sur les solides (DFT périodique). Elle est aussi depuis une dizaines d'années, beaucoup appliquée aux molécules.
Trois méthodes principales existent.
\\
{\GOWO}, {\evGW} et {\qsGW} \tabref{tab:methList}.
Nous nous concentrerons ici sur {\GOWO} et {\evGW}.
\begin{table}
\caption{\label{tab:methList} Liste des méthodes {\GW}}
\begin{ruledtabular}
   \begin{tabular}{ccc}
		 \GOWO & GW à l'itération 0 & Un seul cycle {\GW}\\
         \evGW & eigenvalues GW & auto-cohérent sur les énergies\\
 		 \qsGW & quasiparticle self-consistent GW & auto-cohérent sur les énergies et les orbitales
   \end{tabular}
 \end{ruledtabular}
\end{table} \\
L'équation quasi particule relie l'énergie de l'orbitale  en fonction de $\osigc$. \\
Cas de {\GOWO} (perturbatif)
\begin{subequations}
\label{eq:QP}
\begin{align}
\label{eq:QPO}
\epsilon_l^\text{QP}&=\epsilon_l^\text{HF}+Z_l \qty(\epsilon_l^\text{HF})\sigc{l}\qty(\epsilon_l^\text{HF})
\intertext{Cas auto-cohérent ({\evGW} ou {\qsGW})
dans ce cas il n'est plus nécessaire d'utiliser $Z_l$ car l'auto cohérence corrige ce problème}
\label{eq:QPSC}
\epsilon_l^\text{QP}&=\epsilon_l^\text{HF}+\sigc{l}\qty(\epsilon_l^{QP-1})
\end{align}
\end{subequations} où $\epsilon_l^{QP}$ est l'énérgie quasiparticule, c'est à dire l'énergie obtenue par la méthode $GW$ pour l'orbitale l la ou $\epsilon_l^{QP-1}$ et celle de l'itération précédente.
\colorbox{yellow}{La partie réelle de $\sigc{l}(\omega)$} étant obtenue par  l'intégrale analytique sur les éléments diagonaux de $\osigc$: \cite{Van_S_2015}
\begin{equation}
	\label{eq:SigR}
	\Re{\ev{\osigc(\omega)}{l}} 
	= \sum_{m} \qty[ 
    \underbrace{\sum_{i} 
    \frac{\abs{(il|\vb{\rho_m})}^2 }
    {\omega - \epsilon_i + \Omega_m}}_{\text{Somme sur les occupées}}
    + \underbrace{\sum_{a} 
    \frac{\abs{(al|\vb{\rho_m})}^2}
    {\omega - \epsilon_a - \Omega_m}
    }_{\text{Somme sur les virtuelles}} ]
\end{equation}
avec $\Omega_m$ l'énergie d'excitation et $\vb{\rho_m}$ le vecteur d'excitation.\\
Or, l'équation \eqref{eq:SigR} a plusieurs solutions. Pour savoir quelle est la solution à garder (dite solution quasiparticule) nous définissons
\begin{equation}
\label{eq:Z}
Z_l\qty(\omega)=\frac{1}{
1-\pderivative{\omega}
\ev{\osigc(\omega)}{l}
}
\end{equation}
avec $0 \le Z_l \le 1$.
Nous cherchons ainsi parmi les solutions $s$ de l’équation {\eqref{eq:QP}}la valeur de $Z_l$ la plus grande.
Toutes les autres solutions sont dites satellites.
On a aussi pour $\nsol$ le nombre de solution $\omega_{l,s}$ de l'équation quasiparticule avec $\omega_{1,s}=\epsilon_l^\text{QP}$
et $Z_{l,s}=Z_l(\omega_{l,s})$
\begin{subequations}
\begin{align}
&\sum_{s=1}^{\nsol} Z_{l,s}=1
\\
&\sum_{s=1}^{\nsol} Z_{l,s}\omega_{l,s}=\epsilon_l^\text{HF}
\end{align}
\end{subequations}
\subsection{Problématique}
\label{PB}
L'équation \eqref{eq:SigR} est une combinaison de deux fonctions de la forme $\frac{y}{x}$ avec $x=(\epsilon_n - \epsilon_k \pm \Omega_m)$.
Cette fonction n'est donc pas définie pour $(\epsilon_n - \epsilon_k\pm \Omega_m)=0$
Ce qui entraîne le problème de la discontinuité dans la surface d'énergie potentielle.
\section{PyDuck}
{\PyDuck} est un script utilisant {\quack} pour réaliser des calculs de chimie quantique.
Il se compose d'un script {\python} et dans le même dossier 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.
Ce fichier 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 additionnels pour savoir sur quel graphique l'afficher ou quels orbitales prendre 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 "SinglePoint" où nous pouvons faire des tracés en fonction de $\omega$ ou "Scan" ou les tracés se font en fonction du rayon $R$ pour les diatomiques ou de $Z$ le numéro atomique dans le cas des monoatomique. Les molécules a plus de deux atomes ne sont pour le moment pas supporté dans ce mode.
\begin{comment}
\begin{figure}
        \begin{algorithmic}
                \State Load JSON options
                \if{Option[Mode]==Modes.Scan}
                \FOR R:=Option[Sca] \TO   \STEP OPtion[Scan] \DO
     |expt|(2,i); \\ |newline|() \OD
                \ElseIF{Option[Mode]==Mode.SinglePouint]
                
                \EndIF
        \end{algorithmic}
        \caption{
                \label{fig:algo_PyDuck}
                Pseudo-code for the \PyDuck script.
        }
\end{figure}
\end{comment}
\section{Résultats}
\subsection{Scan}
En faisant des scans nous pouvons identifier différentes molécules où les discontinuités sont susceptible d'apparaître.
\begin{itemize}
\item \ce{H2}.
\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 G0W0@HF 6-31G sur \ce{H2} les trois premier calculs sont réalisé à la géométrie d’équilibre. Le derner et réalisé à la discontinuité LUMO+1}
\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} les trois premier calculs sont réalisé à la géométrie d'équilibre. Le dernier et réalisé à la discontinuité LUMO+1}
\end{figure}
\begin{figure}
\includegraphics[width=0.49\linewidth]{fig/eqpH2evGWzoom1}
\includegraphics[width=0.49\linewidth]{fig/eqpH2evGWzoom2}
 \caption{
        \label{fig:H2evGWZoom}
        Calculs evGW@HF 6-31G sur \ce{H2} zoomé}
\end{figure}
Dans le premier cas nous avons deux divergence de $\osigc$ alors que le deuxième cas est plus classique avec une seule divergence. \figref{fig:H2evGW} entraînant toutes des discontinuités.
\item \ce{HeH+} On observe ici deux discontinuités une à 3 bohr dans la LUMO+2 et une à 4 Bohr dans la LUMO+1.
\item \ce{LiF} Ici nous avons uniquement une discontinuité dans la LUMO+1 à 2.5 bohr.
\end{itemize}
On remarque ainsi rapidement que plus l'orbitale et haute en énergie plus les discontinuités se présentes tôt. nous remarquons aussi une petite bosse en G0W0 l'osque Sigma diverge \figref{fig:H2G0W0}, là où il y a une discontinuité en evGW \figref{fig:H2evGW} cela correspond a un problème de renormalisation \eqref{eq:QPO} pas parfaitement bien réalisé.
\subsection{Fonctions de $\omega$}
Nous avons calculé les valeurs de la fonction $\sigc{l}(\omega)$ ainsi que $Z_l(\omega)$ pour les molécules précédemment étudié à la distance où se produit la discontinuité nous voyons ainsi $Z_l$ qui diminue et $\sigc{l}$ qui diverge \figref{fig:H2evGW} nous avons l'équation $\omega-\epsilon_l^\text{HF}=\sigc{l}\qty(\omega)$ équivalent à \eqref{eq:QP} symbolisé par l'intersection des courbe de $\osigc$ et de $\omega-\epsilon^\text{HF}$ qui nous montre les solutions de l'équation quasiparticule.
Ainsi lorsqu'il y a une discontinuité nous avons plusieurs solutions toutes avec des Z très petit il y a donc un brutal changement de solution d'où la discontinuités.
Pour la suite de notre étude nous nous sommes basé sur les problèmes rencontré par \citet{Van_S_2015}
En effet dans cet article ils font remarquer que sur certaines molécules la discontinuité est présente dans la HOMO au niveau de la géométrie d'équilibre. Ce qui pose un grave problème car dans ce cas nous ne pouvons pas nous fier à ce résultat non physique.

Une première chose que l'on remarque est que le gap HOMO/LUMO DFT est bien inférieur au gap HOMO/LUMO Hartree-Fock. \tabref{tab:gap}
Or, nous avons une grande stabilité de la méthode GW à $30$ fois la valeur du gap ($E_\text{g}$) autour de $0$. Ce qui fait un grand intervalle de valeur stable autour de zéro. Mais nous savons que dans la DFT le gap est très petit, donc si la solution de l'équation quasiparticule se trouve en dehors de la limite des $30 \times E_\text{g}$ les problèmes de discontinuités peuvent apparaître.
\begin{figure}
\includegraphics[width=0.49\linewidth]{fig/SigCO3evGW}
\includegraphics[width=0.49\linewidth]{fig/ZO3evGW}
\caption{
   \label{fig:H2G0W0}
        Calculs evGW@HF 6-31G sur \ce{O3}. réalisé à la géométrie d’équilibre.}
\end{figure}
\section{Solution}
Pour corriger le problème présenté à la section \ref{PB} nous pouvons tenter de modifier la fonction $\sigc{l}$ nous voulons donc changer l'expression située dans l'opérateur somme dans l'équation \eqref{eq:Sig}
Une solution à été testé et semble correcte
les critères pour trouver une bonne fonction $\sigc{l}$ sont :
\begin{itemize}
\item $\sigc{l}(\omega)$ continue et dérivable sur $\numsys{C}$
\item $Z_l\qty(\omega) \in \qty[0,1]$ avec $Z_l\qty(\omega)$ définie à l'équation \eqref{eq:Z}
\end{itemize}
Nous avons testé l'équation suivante
\begin{equation}
\newcommand{\x}{\omega-\epsilon_l\pm\Omega}
\newcommand{\y}{\abs{(kl|\rho_m)}}
\newcommand{\re}{\frac{\x}{\y}}
\y^2\frac{\x}{\qty(\x)^2+\eta^2}
\end{equation}
Dans ce cas, nous pouvons remarquer que sur certaines valeurs de $\omega$ $Z_l(\omega)$ peut être supérieur à $1$
en effet d'après l'équation \eqref{eq:Z} on a
\begin{equation}
Z_l<1\equiv \sigc{l} \textnormal{décroissante}
\end{equation}
\section{Conclusion}
Les méthodes {\GW} posent un réel problème du fait que la divergence de la self-energy induit des discontinuités dans les surface d'énergie potentielles nous avons pu montrer que plus l'orbitale est haute plus les discontinuités arrive tôt. L'autre chose importante est que la méthode est stable au alentour de 30 fois le gap HOMO/LUMO et comme le gap HOMO/LUMO et sous estimé en DFT alors qui est surestimé en Hartree-Fock ce qui fais de Hartree-Fock une meilleur méthode pour éviter ces discontinuités
\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 Bouammali pour sa gentillesse et son aide ponctuelle.
\bibliographystyle{plainnat}
\bibliography{Rapport}
\appendix
\end{document}