Étude ab initio des glaces d’ammoniac fluoré et hydraté sous conditions thermodynamiques extrêmes
Phonons et spectroscopie vibrationnelle
Nous avons présenté le cadre de la théorie de la fonctionnelle de la densité qui permet de déterminer la structure électronique du système pour une configuration donnée des noyaux. En plus des calculs de structures électroniques, nous nous sommes concentrés sur l’étude des vibrations du réseau cristallin et la simulation de spectres infrarouge et Raman. Pour accéder aux propriétés vibrationnelles des glaces d’ammoniac fluoré et hydraté sous pression, nous avons effectué 16 Chapitre 1: Méthodes de modélisation des calculs de spectres afin de guider les expériences sous enclumes de diamant et pour comparer nos spectres aux données expérimentales. Dans la section 1.3.1, nous allons présenter le cadre général du calcul des modes propres de vibration pour un nombre Nn de noyaux, puis dans la section 1.3.2, nous appliquerons le calcul des modes propres de vibration au cas des cristaux. Enfin dans la section 1.3.3, nous montrerons comment, à partir d’un calcul statique (0 K), il est possible de simuler un spectre infrarouge où Raman.
Décomposition en modes propres
A température nulle, au sein d’un cristal ou d’une molécule, tous les atomes sont immobiles à leur position d’équilibre, {R0}. À température non nulle, les atomes vibrent autour de cette position d’équilibre. En réalisant un développement de Taylor limité au terme du second ordre de l’énergie totale E({R}) en terme des déplacements atomiques U α I = Rα I − R (0)α I , nous obtenons : E({R}) = E({R})0 + X I,α » ∂E({R}) ∂Rα I # U α I + 1 2 X I,α,J,β » ∂ 2E({R}) ∂Rα I ∂Rβ J # U α I U β J (1.54) où les lettres α, β désignent les coordonnées cartésiennes (α = {x, y, z}), U α I est la norme du vecteur déplacement de l’atome I par rapport à sa position d’équilibre et E({R})0 est une constante. À l’équilibre le terme d’ordre 1 est nul, ainsi l’Éq. (1.54) devient : E({R}) = E({R})0 + 1 2 X I,α,J,β » ∂ 2E({R}) ∂Rα I ∂Rβ J # U α I U β J (1.55) Dans l’Éq. (1.55), limiter le développement de l’énergie au terme du second ordre constitue l’approximation harmonique. Le terme d’ordre 2 noté C αβ IJ est nommé matrice des constantes de force interatomiques : C αβ IJ = » ∂ 2E({R}) ∂Rα I ∂Rβ J # (1.56) Cette matrice est de taille 3Nn x 3Nn et elle est symétrique, à savoir que : C αβ IJ = C βα JI (1.57) La somme des éléments de matrice sur une ligne quelconque vaut zéro : X J,β C αβ IJ = 0 ∀ I, α (1.58) Cette propriété dérive du fait que, dans un système isolé non soumis à une contrainte extérieure, une translation globale et rigide de la position des atomes ne fait pas varier l’énergie interne du système en question. Si on assimile les noyaux à des particules classiques, leur mouvement est régi par le principe fondamental de la dynamique : MIU¨α I = F α I = − X J,β C αβ IJ U β J (1.59) Par conséquent, il faut résoudre un système de 3Nn équations différentielles linéaires du second ordre couplées entre elles. On cherche des solutions sous la forme d’ondes sinusoïdales : U α I (t) = 1 √ MI u α I e iωt (1.60) Une solution de ce genre est un mode propre : tous les atomes oscillent avec la pulsation ω, à part le facteur 1/ √ MI , u α I représente l’amplitude d’oscillation. Par substitution de l’Éq. (1.60) dans (1.59) on obtient : X J,β 1 √ MIMJ C αβ IJ − ω 2 δIJ δαβ UJ,β = 0 (1.61) D’après le théorème de Hellmann-Feynman, la matrice C αβ IJ peut également s’exprimer comme : C αβ IJ = ∂ 2E({R}) ∂Rα I ∂Rβ J = − ∂F α I ∂Rβ J (1.62) Ainsi, nous sommes partis de Nn équations couplées et nous avons décomposé le problème en Nn oscillateurs indépendants (chaque oscillateur possède sa propre pulsation). L’Éq. (1.61) admet comme solution l’ensemble des coefficients {U β J } et pour déterminer les modes propres de vibration il faut résoudre l’équation séculaire correspondante 4 : det|C αβ IJ − ω 2 Id| = 0 (1.63) où Id est la matrice identité. Comme écrit dans l’article de Baroni et al. [31], il existe plusieurs méthodes permettant de calculer la matrice des constantes de force interatomiques. Or, si on s’intéresse à un solide cristallin, ce dernier possède des symétries ponctuelles (rotation, plan miroirs, …, etc) qui peuvent être utilisées pour réduire le nombre d’éléments indépendants de la matrice C αβ IJ . Dans la section suivante, nous montrons comment utiliser la propriété de périodicité du réseau cristallin pour résoudre l’Éq. (1.61) et trouver l’ensemble des fréquences propres du système. Avec l’utilisation des conditions aux limites périodiques, on ne considère dans la suite que les vibrations au sein d’un cristal parfait, qui est complètement défini par la donnée de sa structure cristallographique (paramètres de maille et positions atomiques).
Modes propres de vibration dans un cristal
Dans cette section, on souligne les spécificités des modes propres de vibration dans un cristal en utilisant la symétrie de translation du réseau. Dans un cristal parfait, le vecteur position du I eme ` noyau peut être développé comme : RI = Rl + τ s = l1a1 + l2a2 + l3a3 + τ s (1.64) avec Rl un vecteur qui repère la position de la l eme ` maille au sein du réseau de Bravais et τ s est le vecteur repérant la position du s eme ` noyau dans la maille primitive. Le vecteur Rl peut être exprimé en fonction des trois vecteurs de base de la maille primitive a1, a2, a3 où les coefficients l1, l2, l3 prennent des valeurs entières. Les déplacements atomiques seront notés U α s tels que : U α s = R α s − R (0)α s (1.65) Avec les nouvelles notations utilisées, la matrice des constantes de force interatomiques devient : C αβ st (Rl , Rm) = » ∂ 2E({R}) ∂Rα s ∂Rβ t # (1.66) 4. En pratique, pour connaître les fréquences propres, on calcule les éléments de matrice de C αβ IJ . Soit il existe une forme analytique simple pour E({R}) et le calcul est possible. Dans le cas contraire, on utilise la méthode des différences finies pour déterminer les forces qui agissent sur chaque atome en supposant que ces derniers effectuent des petits déplacements autour de leur position d’équilibre. Enfin, on diagonalise la matrice des constantes de force pour déterminer les fréquences propres et les vecteurs propres associés. 18 Chapitre 1: Méthodes de modélisation Toutefois, cette matrice a une propriété de symétrie remarquable liée à la symétrie de translation du réseau 5 : elle est invariante par translation de réseau (Rl → Rl − Rm). Ainsi, la matrice des constantes de force ne dépend que de la différence Rl − Rm : C αβ st (Rl , Rm) = C αβ st (Rl − Rm, Rm − Rm) = C αβ st (Rl − Rm) (1.67) Par l’utilisation de cette propriété nous réduisons le nombre de degrés de liberté de 3Nn à 3Nn,c, avec Nn,c le nombre d’atomes dans la maille primitive. Par analogie avec le cas d’une fonction d’onde mono-électronique dans un potentiel périodique, nous cherchons ici les solutions sous la forme d’ondes planes : U α s = 1 √ Ms u α s e i(k.Rl−ωt) (1.68) où u α s est l’amplitude de vibration qui est identique pour tous les atomes équivalents par symétrie et ω est la fréquence de vibration qui reste à déterminer. Si on utilise l’Éq. (1.68) dans l’équation du mouvement (1.59) on obtient : ω 2u α s e −iωt = − X t,β 1 √ MsMt D αβ st (k)u β t e −iωt (1.69) avec : D αβ st (k) = X l e ik.RlC αβ st (Rl) (1.70) D αβ st (k) est la matrice dynamique du système et elle est définie comme la tranformée de Fourier de la matrice des constantes de force interatomiques. C’est une matrice de taille 3Nn,c x 3Nn,c ; nous avons donc réduit le nombre de degrés de liberté de vibration au nombre d’atomes dans la maille primitive. La matrice D αβ st (k) est une matrice complexe hermitique, à savoir que : [D αβ st (k)]∗ = D αβ st (−k). Par analogie avec l’Éq. (1.61), la nouvelle équation séculaire s’écrit comme : X t,β 1 √ MsMt D αβ st (k) − ω 2 (k)δstδαβ U β t = 0 (1.71) Cette équation admet une solution non triviale pour les fréquences ω(k) et les amplitudes de vibration u β v si : det|D αβ st − ω 2 Id| = 0 (1.72) La résolution de l’Éq. (1.72) fournit, à chaque point k appartenant à la 1ZB, les fréquences ω (fréquences propres) et les amplitudes de vibration correspondantes (vecteurs propres). Au total, il y a 3Nn,c fréquences propres de vibration à calculer. En fonction du groupe d’espace considéré, certaines fréquences propres peuvent être dégénérées ce qui permet de diminuer le nombre de fréquences à déterminer et par conséquent le temps de calcul. Les modes de vibration dans les cristaux sont appelés phonons. Il s’agit de modes cohérents de propagation d’ondes de vibration avec un temps de vie infini dans l’approximation harmonique. La relation ω = ω(k) est la relation de dispersion. Le calcul de ω(k) comporte deux étapes : construction de la matrice dynamique D αβ st puis résolution de l’équation séculaire. On montre que pour chaque point k appartenant à la 1ZB, il y a plusieurs modes de propagation différents qui se propagent avec ce même vecteur k. On montre que pour trois de ces modes tels que kkk soit petit, la relation de dispersion prend la forme suivante : 5. Les forces interatomiques F α I sont également invariantes par translation du réseau cristallin. ω1(k) ≈ v1kkk, ω2(k) ≈ v2kkk, ω3(k) ≈ v3kkk (1.73) avec v la vitesse de propagation du son. Ces trois modes de propagation 6 sont appelés modes acoustiques et ils sont responsables de la propagation du son dans le solide. Au centre de la 1ZB, c’est à dire au point Γ, leur fréquence est nulle. En plus des modes acoustiques, il existe les modes optiques de plus hautes fréquences. Dans un mode acoustique, les noyaux au sein de la maille se déplacent en phase, la dynamique est dominée par l’interaction entre les mailles du cristal. Dans un mode optique, les noyaux d’une même maille élémentaire vibrent en opposition de phase les uns par rapport aux autres. Ces derniers sont plus généralement associés aux modes de vibration d’élongation (stretching en anglais), de flexion (bending), de balancement (rocking), d’agitation (wagging) et de torsion (twisting) propres aux molécules et qu’il est possible de détecter par spectroscopie Raman ou infrarouge. Dans la section suivante, nous présentons le calcul de spectres infrarouge et Raman.
Calculs de spectres infrarouge et Raman
Pour le calcul des phonons, nous avons utilisé le code PHONON de Quantum-Espresso [30]. Dans cette thèse, les calculs des fréquences propres de vibration ont été effectués sur des cristaux sous pression. Par conséquent, la maille élémentaire et les positions atomiques auront été préalablement optimisées à la pression souhaitée en utilisant le code PW. Les paramètres de maille et les coordonnées atomiques optimisées sont ensuite utilisés pour réaliser un calcul des phonons. 7 En pratique, les modes propres de vibration sont calculés au point Γ, qui est l’origine de la 1ZB. Par conséquent, la matrice dynamique Dαβ sv (k = 0) du système doit être calculée en ce point. Au sein du code que nous avons utilisé, les modes propres de vibration sont calculés à température nulle grâce à la théorie perturbative de la fonctionnelle de la densité (ou Density Functional Perturbation Theory, DFPT en anglais) [31–33]. La transformée de Fourier de la matrice des constantes de force permet d’accéder à la matrice dynamique qui fait l’objet d’une diagonalisation dans l’espace réciproque pour obtenir les vecteurs propres, U α s (k = 0) et les fréquences propres, ω(k = 0). Cette opération est réalisée en utilisant l’exécutable ph.x du code PHONON. Maintenant que les fréquences propres du système sont connues, il est utile de calculer l’intensité infrarouge ou Raman associée à chaque phonon. Commençons par donner l’expression de l’intensité infrarouge que nous noterons IIR. Pour le ν eme ` mode, IIR s’exprime comme : IIR(ν) = X α X s,β Z ∗αβ s U β s (ν) 2 (1.74) où Z ∗αβ s correspond au tenseur des charges effectives de Born [34] de l’atome s. Plus généralement, la charge effective de Born se définit comme la variation de la polarisation P αβ s dans la direction α causée par le déplacement de l’atome s dans la direction β. La DFPT permet de calculer le tenseur des charges effectives [35, 36]. Pour le calcul de l’intensité Raman, revenons sur les mouvements des électrons et des noyaux au sein du cristal. Deux types de mouvement sont possibles : le mouvement des noyaux par rapport aux autres noyaux à la fréquence propre ων et celui des électrons par rapport aux noyaux à la fréquence propre ωe. Soumis à un champ électrique E d’une onde lumineuse de fréquence ωi , seuls les électrons sont suffisament légers pour répondre à la fréquence excitatrice et vibrent alors à la même fréquence de l’onde incidente. Il y a création d’un moment dipolaire électrique P oscillant tel que, pour des champs faibles : 6. Dans les cas simples, il y a un mode longitudinale et deux modes transverses. 7. Il est évident que l’énergie de coupure choisie pour le calcul d’optimisation de la structure cristallographique doit être la même que celle utilisée pour le calcul des phonons. 20 Chapitre 1: Méthodes de modélisation P = T αβ r E (1.75) où T αβ r est le tenseur Raman également appelé tenseur de polarisabilité. L’intensité de diffusion [37,38] du ν eme ` mode s’exprime comme : IR(ν) ∝ (ωi − ων) 4 ων T αβ r (ν) (1.76) Dans l’Éq. (1.76) α et β sont respectivement les directions de polarisation de l’onde incidente et de l’onde diffusée, et ων est la fréquence de vibration du ν eme ` mode. Les symétries cristallines imposent la forme du tenseur Raman (ces tenseurs sont répertoriés dans la référence [39]) et la symétrie des modes de vibration associés. Suivant les composantes nulles de T αβ r , il existera des modes propres de vibration qui seront inactifs en spectroscopie Raman. Le calcul de l’activité infrarouge et Raman des modes propres de vibration est permis grâce à l’exécutable dynmat.x. Ce dernier calcul requiert la connaissance des fréquences propres, du tenseur des charges effectives et du tenseur Raman (ph.x). On obtient alors un fichier de sortie contenant la valeur des fréquences propres avec l’intensité infrarouge/Raman. En fonction du système étudié, certains modes auront une activité infrarouge nulle mais ils seront actifs en Raman et inversement.Comme le calcul est réalisé à température nulle à partir d’une structure dont les paramètres de maille et les positions atomiques ont été optimisés, la largeur des bandes vibrationnelles obtenues n’a aucune signification physique. Les spectres calculés à 0 K dans cette thèse seront donc présentés sous la forme d’un histogramme.
Remerciements |
