Etude par DFT de photocatalyseurs pour des
applications en photodissociation de l’eau

Modélisation de systèmes périodiques à 2 et 3 dimensions

Détermination de la fonction d’onde d’un système périodique Les systèmes étendus que nous allons étudier (semi-conducteurs, électrodes) sont des solides cristallins et présentent une périodicité à l’échelle atomique. Plutôt que de travailler sur un très grand nombre d’atomes, on pourra utiliser cette symétrie translationnelle en définissant des conditions aux limites périodiques (PBC), le nombre d’atomes à étudier sera alors celui de la plus petite unité de répétition appelée maille élémentaire. Le Hamiltonien aura la périodicité du système, et d’après le théorème de Bloch, les fonctions d’onde mono-électroniques s’écrivent sous la forme de fonctions de Bloch i : ψk(r + R) = ψk(r)e ik·R , ou encore : (II.17) ψk(r) = u(r)e ik·r (II.18) où u est une fonction périodique et R n’importe quel vecteur laissant l’Hamiltonien invariant. Le vecteur k est appelé vecteur d’onde de Bloch. L’espace des vecteurs k est également périodique et constitue le réseau réciproque. Comme vu dans le chapitre précédent, dans un solide les orbitales moléculaires laissent la place à une continuité d’orbitales qui forment des bandes d’énergie. Le vecteur k est un descripteur de choix pour ces orbitales. En effet il suffira de balayer la première période de l’espace réciproque (appelée première zone de Brillouin ou zone de Brillouin) pour décrire toutes les fonctions d’onde du système. Le nombre de points k dans la première zone de Brillouin n’en est pas moins infini, mais on pourra l’échantillonner en un nombre fini de points appelés points k. En effet, les propriétés du système comme l’énergie totale convergent lorsque la précision de l’échantillonnage augmente, on pourra donc choisir le nombre de points k utilisés à partir d’un critère de convergence. Monkhorst et Pack ont proposé une méthode d’échantillonnage homogène de cette zone en points k, 14 qui sera celle utilisée dans cette thèse. Dans certains cas, pour tenir compte d’une substitution d’atome partielle ou encore ajouter un adsorbat sur une surface qui ne ressente pas l’interaction de l’adsorbat de la maille voisine, on pourra être amené à démultiplier la taille de la maille élémentaire, pour i. On peut voir cette équation comme une conséquence de la périodicité de la densité électronique, et donc du module de la fonction d’onde 36 Chapitre II Méthodologie former une supermaille. Cette augmentation de taille s’accompagne de la réduction du paramètre de maille du réseau réciproque, on pourra donc se permettre une diminution du nombre de points k inversement proportionnelle à l’augmentation de la taille de la maille. Maintenant que nous avons défini les fonctions d’onde mono-électroniques à étudier, nous verrons comment elles peuvent être calculées. On distingue deux approches : — L’approche localisée, par décomposition en fonctions centrées sur les atomes constitutifs. — L’approche délocalisée, par décomposition en une série de fonctions d’ondes planes stationnaires. Fonctions de base localisées Cette approche est dérivée de l’approche moléculaire dans le cas d’un nombre fini d’atomes. Dans une molécule, on peut décomposer les fonctions d’onde mono-électroniques comme une combinaison linéaire d’orbitales atomiques (méthode LCAO) : ψi = X j cijχj (II.19) Comme il existe une infinité d’orbitales atomiques par atome, cette méthode bien qu’exacte nécessiterait un temps infini. Une approximation supplémentaire sera alors de tronquer la base des coefficients cij (approximation LCAO). Il reste finalement à déterminer la forme analytique des fonctions χj . On peut s’inspirer des solutions exactes de l’atome d’hydrogène et utiliser des fonctions de Slater (appelées STO), mais la manipulation de ces fonctions est coûteuse en temps de calcul, on préférera les approximer par des fonctions gaussiennes (appelées GTO) qui divisent ce temps d’un facteur 4-5. Pour un atome, une fonction gaussienne est une bonne approximation d’une fonction de Slater loin du noyau mais pas à son voisinage (Figure II.2). Pour remédier à cela on pourra réaliser des combinaisons de gaussiennes dites primitives pour former des gaussiennes dites contractées. On nommera par exemple STO-3G une orbitale de Slater approchée par 3 orbitales gaussiennes. L’approche par des bases localisées consistera donc, une fois l’allure des orbitales atomiques définie, à optimiser les coefficients cij dans les calculs SCF. Dans un solide, les orbitales moléculaires sont remplacées par une continuité d’orbitales appelées orbitales cristallines. Pour déterminer les orbitales cristallines, on pourra étendre la méthode LCAO au cas périodique en les écrivant comme des sommes non pas d’orbitales . Figure II.2 – Allure d’une orbitale STO, et fonctions GTO et STO-3G associées. atomiques mais de fonctions ayant la périodicité du cristal : ψOCi (r,k) = atomes X N cN i(k)φN (r, k) (II.20) Chaque fonction φN (r, k) représente un atome donné sur chaque maille du cristal. On peut alors l’écrire comme une fonction de Bloch : φN (r, k) = mailles X R χN (r − R)e ik·R (II.21) où χN est une orbitale atomique que l’on écrira à son tour comme une combinaison de nG fonctions gaussiennes G χN (r − R) = XnG j c 0 j e −αj (r−R) 2 (II.22) Ce type de base a été utilisé dans cette thèse à l’aide du code CRYSTAL.15 Ce code utilise des gaussiennes localisées ce qui présente l’avantage de calculer efficacement l’échange Hartree-Fock permettant d’utiliser des fonctionnelles hybrides comme PBE0 ou HSE06. Fonctions de base délocalisées Pour ce type de fonctions de base, la fonction périodique u(r) n’est pas basée sur une somme de fonctions locales, mais on l’écrira comme une série de Fourier. On pourra donc décomposer les fonctions d’onde mono-électroniques en une série harmonique d’ondes 38 Chapitre II Méthodologie planes stationnaires. u(r) = X k Cke ik·r (II.23) Comme pour l’approximation LCAO, cette somme infinie devra être tronquée. Pour définir un critère de troncation, on utilise le fait que ces fonctions d’ondes ont la même allure que les solutions de l’équation de Shrödinger pour un électron libre en conditions aux limites périodiques : ϕ(r) = Aeik·r (II.24) k = ~ 2k 2 2m (II.25) La taille du jeu de bases est alors définie par la donnée d’une limite énergétique kmax (energy cutoff en anglais) à partir de laquelle les harmoniques ne seront pas considérées. Dans le cas de fonctions d’onde très localisées, il faudra faire appel à un grand nombre d’ondes planes. C’est le cas des électrons de cœur, qui de plus n’interviennent pas dans les mécanismes que nous serons amenés à étudier. Pour gagner du temps de calcul, les électrons de cœur ne seront pas explicitement calculés et seul leur écrantage sera modélisé. Pour ce faire, on aura recours à une méthode appelée PAW (Projector Augmented Wave) 16 , qui remplace le système noyau+électrons de cœur par un pseudo-potentiel. Ce type de base a été utilisé à l’aide du code VASP17–19

Modélisation d’une surface

Une surface est par définition étendue sur 2 dimensions seulement. Pour modéliser une interface (par exemple entre un catalyseur solide et le solvant, ou encore dans une jonction solide/solide), on peut donc conserver 2 des 3 dimensions périodiques d’un solide. On définira alors une direction selon laquelle la surface sera définie, la périodicité selon cette direction sera alors brisée ce qui conduira à une double interface que l’on appellera slab. Nous calculerons l’énergie de surface suivant la formule suivante : Esurf = 1 2 (Etot,slab − n · Ebulk) (II.26) avec Etot,slab l’énergie totale calculée par DFT, Ebulk l’énergie d’une maille dans le système périodique en 3 dimensions, et n l’épaisseur du slab en nombre de mailles. Si le slab est trop fin, l’influence d’une interface sur l’autre se fera ressentir, ce que nous cherchons à éviter car l’objectif est de représenter une unique interface. On pourra 39 Curutchet Antton ENS de Lyon donc choisir l’épaisseur de ce slab en faisant un compromis entre cet effet et le temps de calcul qui augmente avec le nombre de couches d’atomes (puisqu’on n’a plus de périodicité dans cette direction). Si supprimer une périodicité ne pose pas de problème dans une approche localisée, l’approche délocalisée ne le permet pas, par construction. On peut contourner ce problème en définissant le slab par un empilement périodique slab-vide-slab-vide. . . Dans ce cas il faudra, au même titre que l’épaisseur de surface, choisir une épaisseur de vide telle que l’interaction entre deux slabs séparées par du vide soit négligeable.

Inclusion de la solvatation

Une surface exposée au vide n’aura pas les mêmes propriétés qu’une surface plongée dans un solvant. Il faudra donc prendre en compte l’influence du solvant sur la structure électronique. Ceci peut être fait de façon explicite, en incluant des molécules de solvant dans le calcul, ce qui implique des temps de calculs considérables. Une alternative est la solvatation implicite, développée par Tomasi et. al., 20,21 qui consiste à décrire le solvant comme un milieu diélectrique continu polarisable (PCM). L’énergie d’interaction (purement électrostatique) entre le solvant et la surface est incluse dans la procédure SCF. Dans ce travail de thèse seule la solvatation implicite a été utilisée, à l’aide de l’implémentation dans VASP VASPsol.Cette méthode décrit la constante diélectrique comme une fonction continue vérifiant l’équation de Poisson-Boltzmann (voir Figure II.9). II.3 Retour sur le mouvement des noyaux Les deux sections précédentes II.1 et II.2 nous permettent de déterminer la structure électronique de tous les systèmes que nous serons amenés à étudier. Cependant dans cette détermination, les positions des noyaux ont été traitées comme des paramètres constants (du fait de l’approximation de Born-Oppenheimer). Or, les noyaux se réorganisent dans une échelle de temps plus grande en fonction des contraintes extérieures pour minimiser l’énergie du système. Cet effet se doit d’être pris en compte pour décrire précisément les positions ri des atomes i du système (qui est appelé géométrie du système), et caractériser un état stationnaire (minimum local ou encore point-selle). 

L’optimisation de géométrie

Le balayage de toutes les géométries possibles pour déterminer celle d’énergie minimale serait trop coûteux (rappelons que pour chaque géométrie une convergence de cycle SCF 40 Chapitre II Méthodologie devra être faite), on préférera une méthode de convergence de géométrie à partir d’une géométrie d’essai. Tout comme pour la méthode SCF, nous pouvons établir un cycle qui conduit progressivement à la géométrie optimale à partir d’une géométrie d’essai en y associant un critère de convergence énergétique. Cette méthode de convergence permet de déterminer la position des atomes dans leur géométrie optimale, mais il peut être également intéressant d’étudier les vibrations des atomes constituants le système autour de leur position d’équilibre.

Le calcul vibrationnel

L’étude vibrationnelle du système nous servira d’une part à déterminer une composante de la constante diélectrique d’un solide (voir section II.4.2), et d’autre part à étudier la réactivité de la surface d’un catalyseur (voir section II.5.3). L’objectif est de décrire l’allure de la surface d’énergie potentielle autour d’une géométrie d’intérêt. On peut écrire l’énergie sous la forme d’un développement de Taylor : E(r) = E(0) + X i ∂E ∂ri ! ri + 1 2 X j,i rjHj,iri + . . . (II.27) où E(0) est l’énergie de la géométrie choisie, ∂E ∂ri son gradient pour la coordonnée ri , et Hj,i = ∂ 2V ∂rj ri la matrice Hessienne (dérivée seconde de l’énergie par rapport aux coordonnées des atomes). Dans le cas d’un état stationnaire (minimum local d’énergie, état de transition), les gradients ∂E ∂ri s’annulent, et on assignera un modèle harmonique aux vibrations autour de la position d’équilibre. Par diagonalisation de la matrice Hessienne on pourra déterminer les modes normaux de vibration, et ainsi caractériser l’état stationnaire ou encore remonter au spectre IR, ou à la constante diélectrique vibrationnelle.