List of figures

 4Heat transfer as a function of the prefixed average fictitious kinetic energy E kin of the electronic wave functions for solid aluminum at the melting point. The arrow indicates the kinetic energy required for adiabatic motion of the electrons according to eq. [1.72]. Taken from ref. [75] . . . . 35 v 1.5 Conserved total energy E cons as computed from Eq. [1.46] from Car-Parrinello (CPMD) and Born-Oppenheimer (BOMD) molecular dynamics for a system of GeSe 2.1 Schematic of the three elastic phases observed in network glasses as a function of increasing connectivity. In the cases of truly random networks the intermediate phase collapses yielding an elastic phase transition from a floppy to a stressed rigid. Taken from [1]. . . . . . . . . . . . . . . 2.2 Upper curve, static structure factor for 63-atom model. The lower curve and points are experimental results from ref [6]. Figure taken from [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. 3 Faber-Ziman partial structure factors for liquid GeSe 2 : GGA calculations (solid line), LDA calculations (dotted line), and experiment (dots with error bars) Ref. [12] . S LDA GeGe (k), S LDA GeSe (k), and S LDA SeSe (k) have been shifted down by 2, 1, and 1, respectively. From ref [11] . . . . . . . . . . . . .

Contour plots of valence electronic charge density for a selected Se-Ge-Se trimer within liquid

GeSe 2 in the LDA (upper left) and GGA (upper right). A difference plot (GGA -LDA) is also shown (lower). The scales are in atomic units. From Ref. [13] . . . . . . . . . . . . . . . . . . . . 2.5 Calculated neutron structure factor for liquid GeSe 2 , obtained within GGA (solid line) and LDA (dots), compared to experimental data (circles) ref. [14]. For clarity, LDA curve is displaced downward by 0.4. From ref. [15] [12] (T = 1050 K) and [20] (T = 1373 K). Solid line: theory, [9] (T = 1050 K) and [19] (T = 1373 K). From Ref. [19] . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 S tot neutron (k) of l -GeSe 2 (solid lines) obtained by Fourier integration of the calculated PDF within given integration ranges 0-r c . The dotted line corresponds to the S tot neutron (k) directly calculated in k space. From ref. [11] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Partial pair-correlation functions for liquid GeSe 2 : black line: present BLYP results, red line: PW results of Ref. [11] , open circles: experimental results of Ref. [12]. From Ref. [26] [41] ; third-order isothermal Birch-Murnaghan equation of state fit to the experimental compression data (solid curve). From Ref. [40] 2.13 Dominant densification mechanisms in GeSe 2 glass at high pressure. From Ref. [42] . . . . . . .

Variation of the average coordination number of the 1st peak in g(r) glasses as a function of

pressure. The open circles are glassy GeSe 4 data obtained from integrating the plotted function between the two minima either side of the first peak. The filled circles are GeSe 2 data from Ref. [40].

The blue triangles are ambient data from Ref. [44]. Figure from Ref. [43] . . . . . . . . . . . . .

Partial pair correlation functions for glassy GeSe

2 . The experimental results of Ref. [5] (solid curve) are compared with the results of Ref. [1] with ρ a (dashed red curve) and with the new simulation at the experimental density ρ exp (dotted blue curve with stars) and the present relaxed structure results at ρ b (NVT ensemble, dashed-dotted green curve). . . . . . . . . . . . . . . . . 3.2 The Bhatia-Thornton concentration-concentration partial structure factor S cc (k) for glassy GeSe 2 .

Experimental results of Ref. [5] (black solid line), simulation results with ρ a Ref. [1](red dashed curve), simulation results with ρ exp (blue dotted curve and simulation results of the calibrated system at ρ b (green dashed curve). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. 3 The pressure-volume equation of state for GeSe 2 glass under compression where P denotes the pressure, V is the volume at pressure P and V 0 is the volume under ambient conditions. The measured data from Mei 6 The G(r) functions from experiment (solid dark (black) curves) were obtained by Fourier transforming the measured F (k) functions shown in Fig. [3.4] via the application of a Lorch modification function. The broken dark (black) curves show the extent of the unphysical oscillations at r-values smaller than the distance of closest approach between the centers of two atoms. The solid light (green) curves show the Fourier transforms of the simulated functions shown in Fig. [3.4] 7 The G(r) functions from experiment (solid dark (black) curves) were obtained by Fourier transforming the measured F (k) functions shown in Fig. [3.4] The solid light (green) and dashed line (blue) curves show the Fourier transforms of the simulated functions of simulated systems with and without thermal cycle respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 The pressure dependence of the partial pair distribution functions from FPMD simulations. From the top to the bottom: the partial radial distribution function g SeSe (r), g GeSe (r) and g GeGe (r). . . . 3.9 The ρ/ρ O dependence of the mean nearest-neighbor bond distance r and coordination number n. Results were obtained from (i) neutron diffraction using a N Ge N Se 2 sample on either the D4 [(red)] or PEARL [(blue)] diffractometer, or using 70 Ge N Se 2 and 73 Ge 76 Se 2 samples on D4 and averaging the results [(green)]; (ii) x-ray diffraction [9] 3.11 The pressure dependence of the intra-polyhedral Se-Ge-Se and the inter-polyhedral Ge-Se-Ge bond angle distributions. NVT simulations (left) and NPT simulations (right). . . . . . . . . . . . . . . 3.12 The ρ/ρ 0 dependence of the intra-polyhedral Se-Ge-Se ((red) with vertical error bars) and inter- polyhedral Ge-Se-Ge ((black) with vertical error bars) bond angles as estimated from the measured Ge-Se, Ge-Ge and Se-Se distances. The experimental results are compared to the mean values θ SeGeSe (broken (red) curve with symbols) and θ GeSeGe (broken (black) curve with symbols) taken from the FPMD simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 3. 13 The reduced density dependence of the q-parameter distribution as obtained from simulations of GeSe 2 glass. At each density, the distribution for all Ge atoms is broken down into its contributions from n-fold coordinated Ge atoms (n = 3, 4, 5 or 6). The vertical (red) arrows mark the q value expected for trigonal bipyramidal units. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Pressure dependence of the n-membered rings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 The pressure dependence is given for the fractions of n-fold coordinated Ge and Se atoms, along with the fractions of these n-fold coordinated Ge and Se atoms that contain homopolar bonding defects. The error bars (usually smaller than the symbol size) were calculated according to Ref. [16]. Left [23], red square are simulation data. Green-dashed line represent the third-order Birch-Murnaghan fit to the experimental data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 (a) The FSDP position (FSDPP) as obtained from direct calculation of the total structure factor (green triangle up) and experiments (solid black circles). The principle peak position (PPP) as computed from direct calculation of the total structure factor (magenta triangles down), from Fourier transformed total structure factor (green left triangles), and experiments (blue squares). In addition the PPP height are provided: as computed from Fourier transformed total structure factor (red plus) and experiment (black stars). (b) The pressure evolution of the X-ray structure factor from FPMD direct calculation on reciprocal space (light green lines), from Fourier transform of the partial real space counterpart (blue dashed lines) and X-ray diffraction from Ref. [23] g GeSe (r), g GeGe (r) and total radial correlation function g tot (r). . . . . . . . . . . . . . . . . . . . Experimental results from ref. [5] 4.11 The intra-polyhedral Te -Ge -Te (top panel) and the inter-polyhedral Ge -Te -Ge (bottom panel). T=300K, Simulation BLYP@300K (black lines), simulation BLYP+vdW@300K (red lines) and simulation PBE+vdW@300K (green lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. 12 The fraction of l-fold Ge and Te atoms (l=1, 2, 3, 4, 5 or 6) as obtained from BLYP@300K (black color), BLYP+vdW@300K (red color) and PBE+vdW@300K (green color) models. . . . . . . . .

Top: Ge atoms bonding to homopolar (black circles), corner-sharing (left green triangle) and edge sharing (blue right triangle) configuration. Se atoms belonging to homopolar connections (red square). Bottom: coordination numbers of Ge atoms (red triangles) and Se atoms (green triangle). Mean coordination number from FPMD (blue diamond) and experiments (black stars). . . . . . . 3.22 A representative subset of Ge and Se atoms, in amorphous GeSe

x 5.1 Phase transition from amorphous to crystalline state in a PCM material. A long laser pulse (set operation) with moderate-intensity beam heats the PCM spot above its glass-transition-temperature leading to a high mobility regime that converts to the crystalline energetically favorable phase. A short laser pulse (reset operation) with high-intensity beam heats the crystalline spot above its melting point. By cooling via the contact with the surrounding amorphous, this crystalline spot converts to the amorphous state. The crystalline state is an ordered state characterized by a low resistance compared to the amorphous state, thus a measure of the resistance allows to read the stored information, whether binary "0" or "1". Figure modified from the original figure 1 in [9]. . 5.2 Similarity of some PCMs crystal structures based on cubic structure. refs [21,22] . . . . . . . .

Fragments of the local structure of GST around Ge atoms in the crystalline (left) and amorphous

(right) states. Ref [24] Ring structure transformation in the phase changes crystal-liquid-amorphous. Ref [25]. . . . . .

Ab initio molecular-dynamics model of the amorphous state of the phase-change non-volatile

memory material Ge 2 Sb 2 Te 5 . Ref [28] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Crystallization temperatures T x as a function of composition in the Ga-Sb-Te phase phase diagram.

Ref [42] 4 The partial pair correlation functions: g GaGa (r), g GaSb (r), g GaTe (r), g SbSb (r), g SbTe (r) and g TeTe (r) respectively. Amorphous state (blue lines), liquid state (red line). . . . . . . . . . . . . . . . . . . 12 Snapshot with the main atomistic configurations of glassy Ga 4 Sb 6 Te 3 . . . . . . . . . . . . . . . . 179 6. 13 The partial pair distribution functions relative to model D (blue line), E (red line) and F (green line).180 6. 14 The proportions of l-fold Ga, Sb and Te atoms (l=1, 2, 3, 4 and 5) as obtained for models D, E and F.182 6. 15 molecular dynamics tests and CPU timing tested for a system of GeSe 2 at ambient pressure. Results are displayed for a 10 ps of dynamics on LINUX cluster (Intel processors Xeon L5420) using 48 processors for all tests, the runs were performed using the 3.17 CPMD version. . . . . . . . .

Temperature evolution during the quench of model A ( dashed green line), model B (blue line

First (FPP)

, second (SPP) and third (TPP) peak positions are provided for all models. n αβ corresponds to the partial coordination number obtained by integrating the corresponding peak over the integration range (IR) covering the peak width. Results for previous simulations at ρ a from Ref. [1] and experiment from Ref. [5] compared to new simulation at ρ exp and at the relaxed structure ρ b . 3.2 N Ge-Ge (N Se-Se ) is the percentage of Ge (Se) atoms in Ge-Ge (Se-Se) homopolar bonds, N Ge (ES) is the percentage of Ge atoms forming edge-sharing connections and N Ge (CS) is the percentage of Ge atoms forming corner-sharing connections. Note that in Ref. [6], a molecular dynamics approach was used in conjunction with a reverse Monte Carlo method. )) for Ge-Ge, Ge-Se and Se-Se Faber-Ziman partial structure factors in the total and difference functions (eq. [6.2] to eq. [6.7]). The numerical values take into account the isotopic enrichments of the samples used in the experiments, the coherent scattering length are given by: b 70 Ge I would like to express my deepest appreciation to Mauro Boero not only for giving me valuable feedback, guidance, support and motivation throughout my academic endeavours, but also for his friendship. I owe a large amount of gratitude to Sébastien Le Roux, for being supportive throughout my PhD project, always being there when needed, for his encouragement and understanding. We have had many lively conversations over the past three years and many amazing trips. He has been an excellent friend. Merci Séb! I am very honoured that Prof. Varnek, Prof. Cleri, Dr. Raty, Dr. Boero and Dr. Coasne accepted to read and judge this work and I am grateful for their precious suggestions.

I would like to acknowledge the various collaborators with whom I had the chance to interact. I acknowledge my gratitude to Prof. Phil Salmon and his team, in particular, Kamil Wezka, Keiron Pizzey and Anita Zeidler for the longest fruitful conversations/discussions/suggestions on neutron diffraction experiments.

I would like to express my very sincere gratitude to Prof. Marco Bernasconi for inviting me to Milano and for the collaboration on PCM. I am thankful to Silvia Gabardi for her help and efficiency on performing the SGCPMD calculations.

I am thankful to Christine Tugéne and Dominiquie Quince for their efficient technical help with the IPCMS supercomputer and my own workstation.

In my daily work I have been blessed with a friendly and cheerful office mates, Yu Liu and Anna Demchenko. Thank you for all the time that we spent together. I thank Duval Mbongo Djimbi and Burak Ozdamar for their friendship.

I would like to thank all of my friends in the lab and Cronenbourg campus, who have made the past four years so enjoyable, and without whose friendship and support the PhD would have been a very different experience. This includes, but is not limited to: Saber, Haifa, Sana, Ines, Hanene, Walid, Wafa, Yassine, Anis, Ahmed, Ferdaous, Mohamed, Sarah, Hashim. I would also like to thank my parents and my brother. They were always supporting me and encouraging me with their best wishes.

xvi

Revue de la thèse en langue Française

Revue de la thèse en langue Française

A-Introduction

Parmi la multitude de composés vitreux existant, une catégorie émergente de composés amorphes fait sensation dans la recherche appliquée et fondamentale: les verres de chalcogénures, qui représentent une famille particulièrement intéressante de matériaux vitreux. Leurs propriétés de transparence dans une large gamme du spectre infrarouge assurent leur utilisation dans de nombreux domaines de l'optique et de l'électronique. Leur capacité de stockage de donnée a révolutionné la fabrication des mémoires à accès direct (RAM) et a conduit à une compétition à la recherche de la mémoire universelle. Les verres de chalcogénures sont ainsi des matériaux prometteurs dans de nombreux domaines tels que les fibres optiques, les mémoires non-volatiles ou les capteurs chimiques. Ce travail de thèse s'inscrit dans le cadre des études théoriques ayant pour but l'étude et l'établissement de la structure et des différentes propriétés physico-chimiques des chalcogénures binaires et ternaires sous différentes conditions. Pour ce but, des techniques de modélisation numérique ab-initio ont été employées. En particulier, nous avons utilisé la dynamique moléculaire par premier principes selon l'approche de Car et Parrinello (CPMD) [1].

Ce travail se devise essentiellement en deux parties: La première est consacrée à l'étude des chalcogénures binaires, notamment les verres GeSe 2 et GeSe 4 sous pression et à l'étude des effets des forces de van der Waals et des fonctionnelles d'échange et corrélation DFT sur la structure de l'amorphe GeTe 4 . Dans la deuxième partie, l'intérêt a été porté à l'étude des matériaux à changement de phase pour le stockage de données. Nous avons pu établir la structure à l'échelle atomique d'un nouveau matériau récemment proposé par les expérimentateurs [2], le Ga 4 Sb 6 Te 3 .

xix

B-Méthodes numériques

Dans l'approche CPMD, les orbitales électroniques ψ i sont considérées comme des variables dynamiques classiques qui obéissent à une dynamique électronique fictive. Cette dynamique suit d'une manière auto-cohérente et adiabatique le mouvement des ions, dans le cadre d'une optimisation structurale ou dans celui d'une évolution dynamique. L'énergie d'interaction totale et les forces sont calculées, dans ce contexte, à travers la théorie de la fonctionnelle de la densité (DFT), alors que les équations Newtoniennes du mouvement (pour les degrés de liberté électroniques et pour les ions) sont résolues simultanément. Dans la plupart des cas traités La fonctionnelle BLYP [3,4] a été utilisée pour le traitement de l'échange et corrélation DFT. Les fonctions d'ondes sont développées sur une base d'ondes planes avec un cutoff de Ry. Des conditions aux limites périodiques sont appliquées dans les trois directions et le contrôle de la température est implémenté selon la méthode de Nosé-Hoover [5][6][7].

C-Résultats

Etude de l'amorphe GeSe 2 sous pression.

Notre modèle est constitué de 40 atomes de Ge et 80 atomes de Se. Dans une première étape, notre intérêt s'est porté sur l'amélioration du modèle à pression ambiante. Nous nous sommes interrogés sur les différences existante entre les données expérimentales et les modèles théoriques du GeSe 2 . Ces différences résident essentiellement dans les proportion relatives des motifs structuraux, englobant les tétraèdres GeSe 4 , les liaisons homopolaires, et des liaisons Ge-Se défectueuses trouvées par la diffraction des neutrons. Puisque le modèle simulé à la densité expérimentale (ρ 0.0334 Å -3 ) présente une pression résiduelle de l'ordre de 1 GPa, nous avons procédé à l'élimination de cette pression résiduelle par le calibrage de la densité. Ainsi le système obtenu a une densité légèrement inférieure (ρ 0.0326 Å -3 ) reproduit avec une très bonne précision les données expérimentales comme le montre la figure 1. Dans une deuxième étape, et après avoir amélioré le modèle initial, nous avons soumis le système à des pressions différentes allons jusqu'à 16 GPa, afin de comprendre la variation de la structure avec la pression. Afin de simuler le système sous pression, deux techniques ont été employées. La première consiste à générer une configuration initiale sous pression en donnant à la cellule de simulation des dimensions correspondantes à une densité plus petite que la densité l'initiale. Dans ce cas, la pression est calculée à la fin de simulation et après avoir ramené

Figure 1

Fonctions de corrélation de paires partielles pour l'amorphe GeSe 2 . Voir [3.1] pour plus de détails.

xxi la température du système à T=0 K (simulation NVT). La deuxième consiste à maintenir la pression constante et permettre à la cellule de simulation d'avoir une variation isotropique du volume dans les trois directions (simulation NPT). Chaque système sous pression est soumis à un cycle thermique de durée moyenne de 400ps englobant une phase de chauffage permettant au système de diffuser et d'explorer l'espace des configurations (voir figure 2). Par la suite, une deuxième phase de refroidissement permet l'obtention d'une structure fiable à température ambiante. Les trajectoires collectées à T=300 K ont une longueur typique de 80ps, ce qui conduit à l'obtention de moyennes statistiques à l'équilibre. Couvrant un domaine allant de l'ambiante à 16 GPa, 13 simulations dans l'ensemble NVT et 5 simulations dans l'ensemble NPT ont été réalisées pour permettre de reproduire l'équation d'état mesurée et aller au delà. Nos résultats sont comparés avec la contre partie expérimentale (diffraction de neutrons avec substitution isotopique) sous pression, obtenus récemment au sein de l'équipe du Prof. P.S. Salmon. Afin de valider notre modèle, nous avons calculé et comparé les facteurs de structures totaux à ceux mesurés, en obtenant un très bon accord.

Contrairement à ce qui a été proposé dans la littérature [8], notre étude montre que la xxii densification du GeSe 2 sous l'effet de la pression n'est pas liée à l'augmentation continue du nombre de coordination des différents atomes, le mécanisme de condensation étant ainsi plus complexe. Dans une première étape, et pour des pressions allons jusqu'à 8 GPa, le système résiste à l'effet de la pression et garde sa structure caractérisée par un arrangement tétraédrique des atomes ( 90% des Ge 4 fois coordonnés et 95% des Se deux fois coordonnés) et par un ordre chimique préservé, obéissant à la loi phénoménologique 8-N (n=2.7). La fraction d'atomes de Ge appartenant à des géométries de types edge-sharing et corner-sharing reste autour de 35% et 42% respectivement. En plus, les différentes distributions angulaires gardent la même allure réflétant ainsi les faibles changement topologiques de la structure sous effet de pressions allant jusqu'à 8 GPa.

Figure 3

La variation en fonction de la pression des fractions d'atomes de Ge "Corner-Sharing" et "Edge-sharing", ainsi que la répartition de ce dernier dans des configurations de type Ge xxiii Du côté du mécanisme de la transition, nous avons mis en évidence le rôle particulier joué par les atomes qui forment des liaisons homopolaires. A la première apparition d'atomes sur-coordonnés, l'on observe un maximum de liaisons défectueuses, les liaisons Ge nonhomopolaires redevenant plus nombreuses au-delà de la pression seuil. Par voie de conséquence, ceci nous permet ainsi de conclure que les atomes impliqués dans des liaisons homopolaires sont la force motrice de la transition structurale. Pendant la transition et à partir de 8 GPa, les atomes de Ge de type corner-sharing se transforment progressivement en des atomes de Ge de type egde-sharing. Ceci est accompagné par la formation de plus en plus d'anneaux de 4 et des structures octaédriques, comme il a été montré par l'inspection visuelle de nos modèles ainsi que par le calcul du paramètre d'ordre local. Ce dernier est caractérisé par l'apparition d'un pic autour de la valeur du paramètre d'ordre q=0, typique d'un arrangement octaédrique. L'étude des distributions angulaires montre que les distributions centrées autour de la valeur 109˚, typiques d'un arrangement tétraédrique, changent d'allure sous l'effet de la pression pour se distribuer autour de 90˚. Un deuxième pic apparaît autour de 180˚pour la distribution Ge-Se-Ge en accord avec l'arrangement pseudo-cubique. Les facteurs de structures obtenus par analyse des trajectoires à T=300 K montrent un très bon accord avec les résultats de la diffraction des rayons X publiés par Skinner [9]. Notre étude de la topologie de ce verre sous effet de la pression a montré que le mécanisme de densification est différent de celui du verre GeSe 2 . Les positions des premiers pics dans les fonctions des distributions radiales ne montrent aucun changement important en variant la pression, ceci étant la manifestation que les liaisons interatomiques ne sont pas affectées par la pression. Par contre, les deuxièmes pics subissent des changements importants sous effet de la compression, au niveau de leurs positions et intensités. En particulier, au niveau de la fonction de corrélation g SeSe (r), on assiste à l'apparition d'un nouveau pic à 3.2 Å. Ceci est le résultat du réarrangement au niveau de la deuxième sphère de coordination autour des atomes de Se. La topologie du système est caractérisée par des atomes de Ge 4 fois coordonnés à des Se et des atomes de Se deux fois coordonnés à des Ge, l'ordre chimique étant ainsi préservée car aucune transition n'est détectée. Les distributions angulaires Ge -Se -Ge et Se -Ge -Se ne changent pas sous effet de la pression, sauf pour

Structure de l'amorphe GeTe 4 : Etude des effets de forces de vdW

Les forces de van der Waals (vdW) sont des forces relativement faibles dues à l'interaction à longue portée des électrons non liés, néanmoins elles jouent un rôle important dans le maintien de la cohésion dans plusieurs classes de matériaux (e.g. complexes moléculaires, systèmes étendus, liquides). Nous nous sommes intéressés à l'étude de l'impact des forces de vdW sur les structures des verres à base de Tellure. Dans ce cadre, nous avons aussi réalisé une étude comparative des performances des fonctionnelles "échange et de corrélation" DFT, en particulier la BLYP et la PBE. Pour ce faire, nous avons étudié l'amorphe GeTe 4 en tenant compte des interactions vdW selon le schéma de Grimme [10].

Dans un premier temps, nous avons mis en évidence l'effet des corrections vdW sur la structure du verre pour une fonctionnelle donnée. L'effet le plus remarquable se voit au niveau des distances interatomiques supérieures aux longueurs typiques des premiers voisins et, en particulier, entre les atomes qui ne forment pas de liaisons chimiques (e.g. Ge-Ge et Te-Te).

xxvii Dans ces cas, on assiste, dans les fonctions de corrélation de paires g GeGe (r) et g TeTe (r), à une augmentation de l'intensité des pics associés. Le faible effet sur les corrélations Ge-Te est dû à la faiblesse des liaisons vdW comparées aux liaisons iono-covalentes Ge-Te. Le premier pic de la fonction de corrélation g GeGe (r) simulé avec vdW est beaucoup plus étroit, ce qui induit par conséquence une baisse du nombre de coordination moyen du système comparé aux simulations sans vdW. La prise en compte des forces de dispersion conduit à un meilleur accord avec la fonction de correlation totale et le facteur de structure mesurés expérimentalement. [2]. Le GaST a montré une grande vitesse de cristallisation et une température de cristallisation assez élevée (271˚C) ainsi qu'une basse température de fusion, réunissant ainsi les bonnes caractéristiques d'un PCM.

Dans ce travail de thèse, nous nous sommes focalisés sur la création et la compréhension de la structure de ce matériau dans son état amorphe. Afin de déterminer la densité du système, nous avons simulé (dans le cadre de la dynamique moléculaire par premiers principes et en utilisant le code CPMD) un système de 117 atomes en ajustant la densité à 0.337 Å -3 (6 g/cm 3 ). Après avoir équilibré le système à la température du liquide, la densité a été de nouveau ajustée à la valeur 0.031 Å -3 en minimisant la pression résiduelle. Le liquide a été par la suite refroidi pendant 100 ps à la température ambiante et, en minimisant encore une fois la pression dans la cellule, nous avons fixé la densité de l'état amorphe à 0.0279 Å -3 . Sachant que la vitesse de refroidissement du liquide joue un rôle important dans l'obtention d'une structure, nous avons aussi généré un deuxième amorphe en baissant la température pendant un intervalle plus long de 300ps. xxx La comparaison entre les structures des deux verres a montré qu'une vitesse de refroidissement plus lente est nécessaire pour obtenir un amorphe avec un ensemble de coordination bien défini.

Dans une deuxième étape, et dans le cadre de notre collaboration avec le groupe du Prof. Bernasconi, nous avons pu générer un troisième verre en utilisant la deuxième génération de CPMD (SGCPMD), connue sous le nom de schéma du Kühne [11]. La combinaison des résultats issus des simulations par CPMD et SGCPMD nous a permit d'établir une description précise à l'échelle atomique de la structure de ce verre. Les atomes de Ga, Sb et Te ont des coordinences moyenne de 4, 3.2 et 2.8 atomes respectivement, en très bon accord avec la structure cristalline du GaTe et de l'amorphe GaSb 7 [12]. L'environnement du Ga présente des motifs structuraux de Te et Sb arrangés dans des configurations tétraédriques. Deux géometries peuvent êtres distinguées: la Ga -Sb x Te 4-x et la Ga -(Ga -Sb x Te 1-x ). Le calcul des distributions angulaires montre que les tétraèdres centrés sur Ge présente des angles plus petits que la valeur typique. La très faible proportion des liaisons Sb-Te ouvre la voie à l'existence de deux phases séparées où le système peut être exprimé comme un alliage de (GaTe) 3 -GaSb 6 . Les distributions des anneaux et des cavités dans le matériau sont en faveur de cette séparation de phase.

C-Conclusions

Ce travail de thèse regroupe des études des chalcogénures binaires et ternaires en utilisant des simulations numériques par premiers-principes. Dans la première partie, nous nous sommes intéressés à l'étude des changements des propriétés structurales des verres de la famille Ge x Se 1-x où nous avons modélisé les verres GeSe 2 et le GeSe 4 sous des pressions allant jusqu'à 16 GPa.

Dans le cas du GeSe 2 , nos résultats montrent que la densification de ce verre se déroule en deux étapes séparées par une pression seuile de 8 GPa. On a pu mettre en évidence pour la première fois le rôle des atome homopolaires dans la transition de l'arrangement tétraédrique à l'arrangement pseudo-cubique ainsi que dans le mécanisme de transition.

Dans le cas du GeSe 4 , aucune transformation structurale du réseau n'a été détectée. Nos analyses montrent que ceci est dû à l'existence des chaines de Se qui se plient et se déforment pour remplir les cavités entre les tétraèdres et contribuent ainsi à préserver la structure initiale et l'ordre chimique.

Notre étude du verre GeTe 4 a montrer que la prise en considérations des forces de vdW dans la simulation a un impact important sur la topologie du verre et, en particulier, sur la coordinance des atomes et leur arrangement dans l'espace. Nous avons également montré que le meilleur accord avec l'expérience est obtenu en combinant la fonctionnelle échange et corrélation BLYP avec les corrections de vdW.

La deuxième partie de la thèse a été consacrée à l'étude des matériaux à changement de phase qui peuvent être utilisés comme matériaux de base pour la fabrication des dispositifs de stockage. En particulier, nous avons exploité la méthode traditionnelle de Car-Parrinello et son évolution récente (dénommée "second génération") pour générer et étudier le verre xxxii Ga 4 Sb 6 Te 3 .Nous avons réussi à déterminer la structure à l'échelle atomique de ce verre. Elle consiste à un mélange de motifs structuraux où Ga est majoritairement tétraédrique, Sb coexiste dans des environnements tétraédriques et trigonals alors que Te se trouve principalement dans des conformations trigonal et planaires. Une séparation de phase entre Sb et Te a été suggérée en partant de l'observation que ce verre peut être décrit comme un alliage (GaTe) 3 -GaSb 6 présentant une distribution importante de cavités.

xxxiii

General Introduction

Modelling is a way to explore, understand and explain macroscopic properties of materials by relying on a scheme that is totally devoid of any use of experimental techniques and tools. The atomic scale time evolution of a given network topology can provide information on the precursors of any macroscopic behavior. This level of description allows to track and identify the changes at the atomic scale level at the origin of the measurable/observable properties of a given material. Different levels of atomic-scale modelling exist, depending on the approach used to describe the interatomic interactions. In this thesis, we model disordered materials via first-principles molecular dynamics techniques. By definition, this approach goes beyond classical models, based on empirical interatomic potentials. In the specific case of chalcogenides, a detailed description of the network topology is out of reach of classical models. This limitation is due to the complex ionic-covalent bonding featured by these networks. In amorphous systems (and, in particular, in the binary systems considered in this study), two levels of order could be distinguished: the short range order (SRO) and the intermediate range order (IRO) [1][2][3]. The short range order (2-5 Å) can be defined by three parameters: the number of bonds, the bond lengths and the bond angles. The intermediate range order (IRO) (5-10 Å) is more complex to describe and arises from the occurrence of a regular structural motifs (e.g tetrahedra, octahedra) cross-linked through a specific type of connections. The determination of the structure of disordered systems for which both short and intermediate range occur is a subtle task that goes beyond phenomenological models (such as the continuous random network and the chemically ordered network). In certain cases, the complexity of the network is substantially increased through the existence of unusual bonding types such as homopolar bonds.

In this area of material science, first principles molecular dynamics (FPMD) is a powerful tool to model materials at the atomic scale. Nevertheless, FPMD cannot be considered as a method directly applicable to all systems and to the description of any physical process. From the standpoint of the electronic structure, its fully self-consistent character and the highly non-linear nature of ground state electronic search calls for a huge amount of computational resources. Therefore, for a given system, it is of paramount importance to achieve a good compromise between the affordable resources, the size of the system (intended to be periodic) and the length of the temporal trajectories.

First-principles molecular dynamics, issued from the combination between density functional theory and molecular dynamics, cannot be taken by itself as the expression of an exact theoretical recipe. Density functional theory becomes a predictive tool provided a great deal of attention is devoted to the careful choice of the exchange-correlation (XC) functional employed within the Kohn-Sham Hamiltonian. A large family of recipes are proposed ranging from those based on a conceptual formulation to some others constructed on more empirical grounds. All these issues are addressed in this work to tackle a structural investigation of binary and ternary amorphous systems.

Many disordered chalcogenide systems are well known and have been largely studied, for instance Ge-Se, As-Se, Ge-S, As-S, Ga-Ge-S, Ge-Sb-Te, Ga-Sb-Te among the others. During the last decades interest on the amorphous chalcogenides increased considerably. This interest is not only fundamental but also motivated by strong technological and commercial applications. The most relevant properties of chalcogenide glasses are briefly recalled in the following: Electronic properties; Regarding the semi-conductor character, glass materials have a typical electronic gap of 1 to 3 eV. The metallicity increases respectively for S, Se and Te. Generally, chalcogenide glasses are p-type semi-conductors and feature a reduced electronic mobility. This small conductivity is due to their structural topology [4,5]. In particular, the disordered character of glasses leads to the creation of gap states and band tails. This semi-conductor behavior depends on the external constraints and could be switched to a conductor one through the increase of temperature.

Optical properties;

This class of properties is particularly wide and includes many photoinduced phenomena. Two classes of photo-induced phenomena can be distinguished [6]. The heat-mode, for which the structural changes in the material are induced through the heat generated by the non-radiative recombination of photo-excited carriers. In the Te-based materials, the heat-mode has a spectacular manifestation, since upon laser heating the structure undergo a reversible thermal crystallization [7]. The second mode is the photon-mode, for which illumination lead to a change in the physical properties of the materials. For example, the photo-darkening phenomena, generated through a red shift of the absorption edge during or after light illumination [8]. This shift could be reversed to the high wave lengths through atomic vibration (annealing at the glass transition). The same photon excitations can lead to changes in the density and in the elastic and chemical properties of the sample [9,10]. Furthermore, photo irradiation could lead to non-linear induced properties such as the optical anisotropy. Thermal properties; Chalcogenide glasses are generated by rapidly quenching from the liquid state. The high quenching rate leads to a rapid increase of the system viscosity generating a glassy system. The liquid-glass transition occurs around a specific temperature called the glass temperature T g . Above this temperature, the amorphous system can attain various structural configurations, however, when going below T g , the system is frozen. These properties allowed for chalcogenides to find application in xerography ("dry writing"), photovoltaic, thin film active devices, bio-sensing, optics, telecommunication, photo-sensing, waveguides, IR fibre optics, medical imaging and non-volatile memory devices. The latter application worth of mention has become an exciting field that stimulated the curiosity and the engagement of scientists and engineers, in the search of materials that provide better performances. We are considering here the so-called "phase-change" materials. Ternary Te-based materials are the most used family materials in this area. A non-exhaustive list of references could be found in these books [4,5,7,11,12] and the references therein. The complexity of the amorphous networks can be substantially increased by altering the character of the bonding. In practice, this could be achieved via doping or through the application of an external thermodynamical constrain: temperature or/and pressure. Since 1953, Bridgman and Simon showed the effects of very high pressures on glassy systems. The elevated load leads to a permanent densification phenomena via a significant increase in the density [13]. This permanent densification phenomena has important implications at the fundamental and industrial levels. It allows the design of materials with tunable properties [14]. The different physico-chemical properties of the glass could be substantially changed within the same composition. Many studies investigated the effect of the pressure in the oxide and binary chalcogenide glasses [15][16][17][18][19].

In this context, molecular dynamics simulations come to highlight the related physical phenomena. Given the above context, the present work have two main purposes: First, we attempt to fill the gap between the experimentally proposed densification mechanisms of binary compounds belonging to the Ge x Se 1-x glasses. X-ray diffraction [20], raman spectroscopy [21] and acoustic measurements [22] proposed different, and yet contradictory, descriptions of the densification of GeSe 2 under pressure. Our implication will mainly focus on providing a fully self consistent set of calculations on the behavior of this glass (GeSe 2 ), and more generally, the behavior of Ge x Se 1-x family of glasses under load via the consideration of glassy GeSe 4 . In addition to this (in a separate part of this thesis work), we are aiming at contributing to the ongoing research on improving the modelling of Te-based chalcogenides. First, we contributed to elucidating the effect of the van der Waals (vdW) dispersion forces on the glass structure of GeTe 4 . This is done by considering two different exchange-correlation functionals. An additional purpose of this thesis is to establish an atomic scale description of phase change materials (PCM) by studying a recently proposed PCM, Ga 4 Sb 6 Te 3 . This compound showed interesting properties that could be useful in non-volatile memory manufacturing.

Based on what presented above, this thesis is divided into three main parts, each part containing distinct chapters: Part I: An overview of first principles molecular dynamics methods is provided in chapter 1.

The usefulness of atomic scale modelling is discussed with examples on advantages and pitfalls of this scheme. A brief introduction to the density functional theory (DFT) is proposed in close connection to the simulations that we are aiming at performing. In addition to their practical implementation, molecular dynamic recipes (say CPMD and BOMD) are introduced and their crucial control parameters are discussed.

Part II: State of the art and results on binary chalcogenides simulations and experiments are provided in this part. Chapter 2 summarizes the chronology of the development of the atomic scale modelling that allowed and motivated the present work. In addition, the main results quoted in literature on the pressure-induced changes in GeSe 2 and GeSe 4 systems are presented. In chapter 3 we provide the results on the study of binary chalcogenides, namely the amorphous GeSe 2 and GeSe 4 , under the effect of pressure. First, we show how the experimental equations of state can be reproduced together with the ambient pressure models. Second, we show results on the topological changes affecting the networks under load. A detailed atomic scale description is given through the study of real and reciprocal space properties, with a close comparison to experiments by neutron diffraction with isotopic substitution (NSID). In addition, further insights are provided via the local environment analysis, this allowed us to elucidate the longstanding controversy related to the nature of the density-driven network collapse and the related mechanisms. Chapter 4 bridges part II and part III. GeTe 4 is a particular chalcogenide glass, it finds applications as a phase change material and as a common opto-electronic device. By modelling the structure and taking into account the van der Waals dispersion forces, we highlight the General Introduction importance of these corrections in reproducing the experimental results. Comparisons with models lacking the dispersion forces are also shown. The importance of the XC functionals employed in modelling materials is highlighted through the use of two different ones. To this end, we compare and rationalize the performances of PBE and BLYP functionals.

Part III: Chapter 5 recalls the essence on phase change materials and the physics behind their outstanding behavior. In particular, we give insight into the performances that a universal memory should feature in addition to the reasons behind the high crystallisation speed of PCMs. Previous experimental works and motivations of the present work are presented. In Chapter 6, methodology and sample preparation methods of amorphous Ga 4 Sb 6 Te 3 are detailed. The structure of this PCM has been elucidated by combining the standard (also called "first generation") Car-Parrinello molecular dynamics (CPMD) to the "second generation" one (SGCPMD). The SGCPMD method allows reducing the computational effort via a substantial increase in the dynamical time step; combined with a more rapid search of the converged electronic ground state for each ionic configuration.

Part I

Atomic scale modelling ATOMIC SCALE MODELLING: WHY AND HOW? Summary

In this chapter I introduce the basic notions of molecular dynamics, to be considered as a tool to achieve an atomic-scale description of a condensed matter system. At the beginning, a short survey of classical molecular dynamics will be given, by stressing advantages and pitfalls. In view of the conclusions drawn herein, I shall describe the main ideas of Density Functional Theory (DFT), seen as a method to obtain accurate structural and electronic properties. In addition to a detailed presentation of the Born-Oppenheimer and Car-Parrinello molecular dynamics methods, some hints about the thermodynamical control methods are also given.

Modeling materials by classical molecular dynamics : advantages and pitfalls 1.Why we need atomic scale modeling

Computer simulation is a powerful tool to explore and investigate atomic scale properties of many-body interacting systems and find applications in a wide range of fields, like chemistry, physics and biology. When the targeted degrees of freedom are the atoms (and, whenever needed, the related electronic structure), one can define "computer simulation" as "atomic-scale simulation", making explicit the fact that the description of a given system will take place at the atomic scale. Within this given context, the power of computer simulation rests on two factors. First, one can access the microscopic properties by using rigorous statistical mechanics techniques, thereby ensuring a well founded description, devoid of any "trial and error" empirical character. Second, the prediction of macroscopic properties becomes accessible, in the range and beyond the range defined by pertinent measurements, fixing the boundaries of the experimental counterpart.

Computer simulation methods at the atomic scale level are basically covered by two main techniques: the Monte Carlo (MC) and the Molecular Dynamics (MD) simulations. In the following, we shall focus essentially on MD as the working strategy of this thesis.

Modeling via classical molecular dynamics

The aim of the MD simulations is to reproduce, as a function of time but with no explicit dependence on time of any variable, the evolution of the forces acting among the atoms. Within classical molecular dynamics (CMD), this is possible once suitable analytical expressions for the interaction potential are made available. In classical molecular dynamics, the interaction potential do not depend explicitly on the electronic structure, this meaning that the force field governing the interaction is a function of the interactomic distances (and angles if needed) among the atoms.

Let's consider a system of N interacting particles described by their cartesian positions R I = (X I ,Y I , Z I ), where I = 1, ..., N. These particles interact via a given potential V (R I ) depending on their positions. The forces f I on particle I are simply the gradients of the potential V (R I ).

f I = - ∂V (R I ) ∂ R I (1.1)
If one postulates that the atoms are subjected to move on time according to the classical mechanics laws, the dynamics can be derived from analytical mechanics by writing the Lagrangian depending on the dynamical variable R I (t):

L(R I , ṘI ) = K -V = N ∑ I=1 1 2 M I Ṙ2 I -V (R I ) (1.2)
Where M I is the I th particle mass. The equation of motion could be derived from the following

∂ L ∂ ṘI = ∂ L ∂ R I (1.
3)

The outcome of equation [1.3] is the Newton's equation of motion:

M I RI = -∇V (R I ) = f I (1.4)
Solving this equation provide the time evolution of the atomic trajectories, opening the way to the evaluation of structural and dynamics properties through the use of temporal averages. The a priory unknown form of the model potential makes the resolution not immediate, this calling for the construction of approximate models. As a first approximation (quantitatively suitable for a restricted class of systems only!) one can use for the atom-atom interactions pair-wise interactions. This means that particle 1 interacts with particle 2 through a potential v(r 12 ) where r 12 = |r 1 -r 2 | (and similarly for all pairs). The total potential reads:

V (R I ) = N ∑ I=1 N ∑ J>I v(r IJ ) = 1 2 N ∑ I,J=1 v(r IJ ) (1.5)
Consequently, the force acting on particle I is now a superposition of pair forces and given by:

f I = - ∂V ∂ R I = - N ∑ J =I ∂ v(r IJ ) ∂ R I = N ∑ J =I f IJ (1.6)
In case of a pairwise potential, the force between two atoms depends only on the distance between the atoms. This is not the case when one invokes N-body potentials, for which the force between two given atoms depends also on a certain number of neighbors.

In practice, classical molecular schemes suitable for atomic-scale simulations are based on three steps : First, we need a model of interaction between the atoms. Second, a stable integrator is needed to propagate the particle positions and velocities from time t to t + δt. Notice that the choice of δt is a crucial key to ensure the stability of the algorithm. Third, we need to make a choice for the simulation statistical ensemble (NVT, NPT or NVE)

allowing to obtain macroscopic averages from the time trajectories of microscopic degrees of freedom (the atoms) through averages taken at thermal equilibrium.

The following general scheme illustrates the main steps of classical MD simulation: The strength of the MD method rely on the fact that, once the potential energy function V (r) and the equations of motions have been chosen, the evolution of the system is completely defined. This evolution can be tracked at the smallest relevant time and length scales and all particle correlations are preserved and can be calculated. Effective potentials for classical molecular dynamics can depend on many parameters and could be refined by fitting experimental data or an ab-initio calculation database. Needless to say, the accuracy of the results strongly depends on the reliability of the used model potential.

To provide an example of systems that cannot be treated by using pairwise potentials, one can quote metals [1]. In this case, the use of pairwise additive potential functions leads to incorrect relationships between the components of the bulk modulus [2]. Therefore the simulation of a metal requires a truly N-body potential function.

Why and when classical molecular dynamics works

Classical molecular dynamics (CMD) methods stands as a powerful tool to study and solve a large class of problems on different fields: defects in solids, fracture, surface properties, molecular clusters, polyelectrolytes, soft matter and bio-systems. In all of these issues, the precise description of chemical bonding at a quantitative level is not an issue of paramount importance. As such, other goals are pursued, most of them related to the possibility of dealing with very large systems. The systems of interest can be very different and range from small molecular sys-14 1.1. Modeling materials by classical molecular dynamics : advantages and pitfalls tems (N < 500) to complex glasses and bio-complexes (N ∼ millions). For instance, in the field of material engineering, the ultimate goal is to fabricate and design novel nanoscale systems, although the control and the understanding of their properties is not always an easy task. The need to control the change of some properties such as the elastic moduli, mechanical strength, toughness, durability, call for a scaling down of mechanical engineering concepts. In this spirit, CMD plays an important role. Recent advancement in computational methodology and highly parallel computers allowed researchers to carry out multi-million-atom simulations with typical linear scale range of 50 to 100 nanometers [3][4][5][6][7]. The Silica glass, which shows outstanding mechanical properties, has been the subject of large scale classical MD simulation [6] where more than one million atoms system have been considered. These simulations allowed to access the mechanisms of shear-induced void deformation. High-strength ceramics (B 4C , SiC, Si 3 N 4 ), featuring remarkable properties such as high hardness and strength, have also been studied via CMD in order to understand the microscopic mechanisms of impact damage, which will allow the design of better ceramic compound for coatings and armors. Chemistry is an other challenging field where CMD is highly used to study extremely complex systems at a qualitative level. Molecular dynamics highlights the role of molecular architecture in determining phase morphology and the related physicochemical properties. For the purpose of designing new molecules of interest to industry, many self-assembling systems were extensively studied. For example, the phase transformation from hexagonal to lamellar structure of 800,000 coarse grain particles was achieved by CMD of an aqueous surfactant solution [7].

Failures of classical molecular dynamics: specific cases

Chalcogenides are materials that show a high ability to form networks over a wide range of compositions (network-forming glasses). The complexity of the network makes the understanding of the bonding nature and the local structural arrangement of the different species a real challenge. Over the previous decades, many experiments and theoretical studies tried to unveil the atomistic structure and the physico-chemical properties of liquids and glasses chalcogenides [8,9]. Within the Ge x Se 1-x (0 < x < 1) family of glasses, many systems shows particular features that can not be reproduced by model potentials. For instance, let us consider that the network could be described by a well defined set of interconnected coordination polyhedra. These are taken as responsible of the intermediate range order (IRO), namely a form of order that establishes on lengths of the order of several interatomic distances [10,11]. Among all the concentrations, the x = 0.33 corresponds to GeSe 2 which is considered as a prototypical network-forming glass. This could be understood in terms of interlinked GeSe 4 tetrahedra.

GeSe 2 shows distinctive challenging features such as the existence of homopolar bonds and the presence of the first sharp diffraction peak (FSDP), characteristic of the intermediate range order [10,11]. In order to model consistently the structure of GeSe 2 , Vashishta and coworkers published several results on the GeSe 2 structure [12,13] where the issues of structural arrangement of the network and the origin of the intermediate rang order were addressed. In these works, the structure of GeSe 2 has been discussed via molecular dynamics based on two-body and three-body interparticle potentials where the atomic charges were taken to be +4Z for the Ge and -2Z for the Se with Z = 0.33e. Vashishta's works described correctly the nearest neighbor distances of Ge -Se and Se -Se, and the number of tetrahedra connection, where 32% of Ge atoms were found in edge-sharing conformations (see appendix C for definition). The FSDP was associated to Ge -Ge and Se -Se real space correlations in the range of 4 -8 Å . Interatomic potentials (similar to the one used in this study), were employed also for other materials like silicon dioxide SiO 2 [14].

An important observation stems from the analysis of the papers produced by the Vashishta team. Indeed, whereas the FSDP appears in the total structure factor S(q), it does not appear in the charge-charge structure factor S ZZ (q) regardless of the specific details of the model potential. This motivates those authors to claim that "the absence of the FSDP in S ZZ (q) is independent of the interaction potential used to describe the system, or the particular detail of the connectivity of the system, or the theoretical scheme used to carry out the calculations" [14]. This statement was employed to associate the absence of the FSDP in S ZZ (q) to the absence of charge-density fluctuations in the 4 -8 Å range. It has to be remarked that, within a pointcharge effective potential model, the charge-charge structure factor S ZZ (q) becomes the same as the concentration-concentration structure factor S cc (q). This describes the ordering of the two chemical species where each peak is related to the preference for like or unlike neighbours at a given distance.Therefore, from the point of view of classical molecular dynamics, a general statement on the absence of fluctuations of concentration on the intermediate range distance could be extracted from the results of Vashishta, this leading to profound implication on the physics and chemistry of glasses. Worth of interest is the observation that, at that time, no experimental evidence existed on the shape of the concentration-concentration structure factor for any material. This situation changed completely two years later, showing the limits of the classical molec-17 1.2. Density Functional Theory ular dynamics approach based on point charges. Indeed, a full set of partial structure factors on liquid GeSe 2 obtained by isotopic substitution with neutron diffraction experiment were published by Penfold and Salmon [15]. Despite the agreement between experiment and simulation on many aspects, this work brought certain evidences that the model potential fails to predict many crucial properties. The correct position of the FSDP was found to be shifted with respect to the experimental one. Also the fraction of ES tetrahedra was "seriously underestimated" in the melt state. More importantly, the FSDP in the S cc (q) was present, thereby unsubstantiating the above claim based on the equivalence between S ZZ (q) and S cc (q). Another important point that was discussed by Salmon is the existence of the wrong (homopolar) bonds, in principle not compatible with the notion of charge of equal sign forming a bond. This is a further demonstration of the inadequacy of Vashishta's interaction potential, which is based on formal charges distribution [16] preventing atoms bearing the same charges to stand in nearest-neighbor interactions. As a complement to this description, it should be added that a step toward the improvement of the Vashista's work was attempted by Mauro and coworkers. These authors developed an interatomic potentials for GeSe 2 based on Moller-Plesset perturbation theory and a cluster expansion approach, where two-and three-body interaction potential are developed by fitting continuous functions to the discrete set of ab-initio data [17]. The main improvement is the reproduction of a small fraction (3.1%) of Se -Se homopolar bonds, while the Ge wrong bounds are still out of reach. In summary, disordered network-forming materials (and in particular chalcogenides) show complex features that cannot be described within the model-potential approximation. It appears that, for the determination of the potential energy surface, an explicit account of the electronic structure is needed. This is exactly the reason why the present work is entirely based on an electronic structure approach combining chemical accuracy and access to reasonable system sizes, that is to say, density functional theory (DFT).

Density Functional Theory 1.2.1 DFT in short

Density functional theory provides a powerful tool for electronic structure calculations of atoms, molecules and solids. The initial idea leading to the DFT was conceived by Thomas and Fermi in 1927 [18,19]. Based on Fermi-Dirac statistics, they derived the electronic density of state and the electric field for heavy atoms, by assuming the electron gas as a completely degenerated Fermi gas. In the middle of the sixties, Hohenberg, Kohn and Sham established a rigorous density functional theory of the quantum ground state based on the basis of quantum mechanics [20]. In the following I will sketch the main steps leading to DFT. Hohenberg, Kohn (HK) showed rigorously by mean of two well established theorems that the total energy of a many-electron system in an external potential V ext is uniquely determined by the ground state electronic density n(r):

E[n(r)] = F HK [n(r)] + d 3 rV ext (r)n(r) (1.7)
F HK is the Hohenberg and Kohn universal functional that does not depend on the external potential. In addition, the exact ground state density is given by the density that minimizes the total energy E[n(r)]. As a consequence, henceforth, we could replace the many-particle wave function by the ground state electronic density as a basic quantity to solve the electronic structure problem. One should note that the explicit form of F HK is unknown, this entailing all the quantum complex many-body effects.

In a further step, Kohn and Sham (KS) proposed one year later a way to reduce the problem complexity without any drastic loss of any theoretical accuracy [21]. They proposed that the electronic charge density could be decomposed in a set of one-particle orthonormal orbitals. In other terms, they replaced the original interactive problem by an auxiliary non-interactive one where the instantaneous electrons repulsion is replaced by an effective potential. This holds true provided the two systems have the same electronic density.

n(r) = N ∑ i=1 | ψ i (r) | 2 (1.8)
The proposed version of the functional F HK reads:

F HK = T KS [n(r)] + E H [n(r)] + E xc [n(r)] (1.9)
where T KS [n(r)] is the kinetic energy functional of non interacting mono-electronic orbitals:

T KS [n(r)] = - 1 2 ∑ i ψ * i (r)▽ 2 ψ i (r)dr.
(1.10)

Density Functional Theory

The second term is the electrostatic Hartree term:

E H [n(r)] = 1 2 n(r)n(ŕ) | r -ŕ | drd ŕ. (1.11)
The last term is the exchange and correlation (XC) energy obtained by simplifying the complexity of the unknown part of F HK by subtracting the known contribution. The explicit form of

E xc [n(r)
] is unknown and accounts for all quantum effects.

E xc [n(r)] = T [n(r)] -T KS [n(r)] + E ext [n(r)] -E H [n(r)]. (1.12)
Gathering all the terms, the total energy functional reads:

E KS [n(r)] = T KS [n(r)] + E H [n(r)] + E xc [n(r)] + E ext [n(r)]. (1.13)
where

E ext = -n(r) N ∑ I=1 Z I |r -R I | dr -∑ I<J Z I Z J |R I -R J | (1.14)
The variational principle applied to Eq. [1.13] allow us to get a mono-electronic equation reflecting the non-interacting system in an effective potential:

(- ∇ 2 2 +V e f f (r) H KS )ψ i (r) = ε i ψ i (r), (1.15)
where

V e f f (r) = 1 2 n(r ′ ) | r -r ′ | dr ′ +V xc (r) +V ext (r), (1.16)
the exchange-correlation potential and V ext are defined as:

V xc (r) = δ E xc [n(r)] δ n(r) (1.17) V ext = δ E ext δ n(r) = - N ∑ I=1 Z I |r -R I | (1.18)

Exchange and correlation approximations

As sketched above, DFT is a formally exact theory for the ground state. Nevertheless, the exact analytic form E xc [n(r)] is unknown. Therefore, we need to introduce some approximations to make the theory tractable. Since the beginning of DFT history, many approximations for In this scheme, the ladder bring us from the "earth" (Hartree theory) to the "heaven" (complete final theory) where the exact XC functional does exist. The intermediate crossbars gather functionals according to their complexity level. Moreover, functionals could be divided into:

empirical (derived by fitting known results on atoms or molecules), and non-empirical based on some physical rules. In the following, we will briefly introduce the most common XC functionals over the first crossbars of the ladder.

Local Density Approximation

The simplest approximation for E xc [n(r)] was originally introduced by Kohn and Sham. It consists on assuming that the n(r) changes very slowly. One could assume that n(r) is locally uniform and equivalent to the same quantity as in the homogeneous electron gas ε hom XC :

E XC LDA [n(r)] = R 3 n(r) ε hom XC [n(r)] dr (1.19)
In the case of the homogeneous electron gas, the XC energy is obtained through the parametrization by Perdew and Zunger [23] based on the Quantum Monte Carlo (QMC) simulations results 20 21 1.2. Density Functional Theory by Ceperley and Alder [24]. Since LDA is derived from the homogeneous electron gas, one expects that it should hold only for systems with slowly varying densities. On the other hand, LDA have been extremely successful in predicting molecular geometries [25], vibrational properties [26], single-particle properties [26] and surface diffusion barriers in simple adsorption systems. However, LDA is far from being satisfactory for cases in which the nature of bonding is a blend of ionic and covalent characters. Although LDA provides a good approximation that works well for several systems, it does produce some systematic errors which of course are relevant to try to remove or at least diminish. For instance, LDA underestimates the exchange energy up to 15% and overestimate the correlation energy up to 200% [27,28]. In addition, LDA overestimates the bonding leading to too small lattice parameters and too large bulk moduli [29,30]. LDA has also been used to calculate barriers for activated adsorption, but calculations on H 2 on Cu(111) and Al(110) showed that the calculated barriers are too low compared to the experimental ones.

Generalized gradient approximation

The natural way to refine LDA is to take into account the inhomogeneity of the electronic density. which gives rise to a family of approximations called the Generalized Gradient Approximations (GGA) [23,31]. This could be done in two different ways. The first way is to use the so-called gradient expansion approximations (GEA), where gradient corrections (|∇n(r)|,

|∇n(r)| 2 , |∇ 2 n(r)|
) are added to the LDA. These power-series like gradient expansions do not add much if lower power corrections are included, while higher order corrections are extremely difficult to calculate. The second way consists in defining a new general function that depend on n(r) and ∇n(r), the typical form for the GGA functional is:

E xc [n, ∇n] = f GGA (n(r), ∇n(r))dr (1.20)
Due to the flexibility on the choice of f GGA (n(r), ∇n(r)), and depending on the target properties of interests (structural properties, electronic properties, band structure, etc), many parametrization have been developed for the exchange-correlation functional. Usually the contribution to the exchange and correlation terms are developed separately:

E GGA xc [n, ∇n] = E GGA c [n, ∇n] + E GGA x [n, ∇n] (1.21)
One of the most popular E GGA c correlation function has been derived by Lee, Yang and Parr [32] and known as LYP functional. The remaining E GGA x correlation part is the most challenging.

Within GGA the latest part reads:

E GGA x [n, ∇n] = n(r)ε uni f x (n(r))F GGA x (s)dr (1.22)
Where ε uni f x is the exchange energy per particle of a uniform electron gas. This quantity is well established in the literature [33]. However, the form of

F GGA x is still a matter of debate.
In what follows, s is a dimensionless reduced gradient:

s = |∇n(r)| 2(3π 2 ) 1/3 n(r) 4/3 (1.23)
For the exchange part, the analytical forms developed within the PBE and B88 functional are worth mentioning. Their analytic form are respectively given by:

F PBE x (s) = 1 + k - k 1 + µs 2 /k (1.24) F B88 x (s) = 1 + β x(s) 2 C[1 + 6β x(s) sinh -1 (x(s))]
(1.25) where x(s) = 2(3π 2 ) 1/3 s and C = 3 2 ( 3 4π ) 1/3 . In case of PBE, the functional are built from physical constraints, thus the parameters k and µ are non-empirical. In B88, β parameter is obtained by empirical fitting. Overall, GGA gives good results and enhances the cohesion energy as well as the lattice parameter. Nevertheless, this improvement is not always in the good direction, since sometime the GGA over-correct the LDA results [34][35][36].

Hybrid functionals

Hybrid functionals are the penultimate step in the Jacob's ladder of the exchange correlation functional approximations (Fig. [1.2]). The original idea goes back to Kohn and Sham. They proposed that an account of the correct exchange term based on the Hartree-Fock approximation added to the electron gas correlation energy could lead to a better asymptotic behavior of the effective potential (∝ -1/r) in the region far from the atom [37]. Despite the success of this hybrid methods on studying atoms, it drastically failed for molecule. In addition, this method involve multiconfiguration self-consistent field calculations (MCSCF) [38][39][40][41][42] or generalized valence-bond (GVB) reference states [43]. Another issue related to this scheme has to do with the fact that the exact exchange energy functional is expressed in terms of the Kohn-Sham 22 1.2. Density Functional Theory orbitals rather than the density, which make the calculation very heavy for a periodic system. One possible solution to enhance the efficiency of hybrid functionals is to reformulate the separation between the exchange and the correlation term, by adding only a portion of the exact exchange from Hartree-Fock theory with the exchange from GGA's functional. Many Hartree Fock based schemes have been proposed, the most commonly used are B3LYP [32,[43][44][45], PBE0 [46], TPSS [47] and the HSE family [48][49][50][51][52][53][54]. During the course of this thesis I had the opportunity to take advantage of one specific hybrid functional (the HSE one) in an attempt to improve the performances of DFT. In what follows an illustration of the HSE theoretical framework is provided.

The HSE exchange-correlation functional is due to Heyd, Scuseria and Ernzerhof. It employs an error function screened Coulomb potential to calculate the exchange portion of the HF interaction energy. The spatial decay of the HF exchange interaction speeds up considerably by changing the 1/r Coulomb potential with a screened one that is split, in the HSE case, into short-range (SR) and long-range (LR) components:

1 r = 1 -er f (ωr) r SR + er f (ωr) r LR (1.26)
where ω is the screening parameter that defines the separation range. The proposed recipe of the exchange correlation energy is:

E HSE xc = aE HF,SR x (ω) + (1 -a)E PBE,SR x + E PBE,LR x + E PBE c (1.27)
where E PBE,SR

x and E PBE,LR x are the short and the long range order PBE exchange functional,

E HF,SR
x is the short range HF exchange and E PBE c is the PBE correlation energy. The parameter a = 1/4 is the HF mixing constant obtained from perturbation theory [55]. For a vanishing screening parameter ω, the HSE reduces to the hybrid PBE0 functional while for an infinite value HSE is nothing else than PBE [55]. Thus, a fixed non zero value of the screening parameter leads to a functional that could be read as an interpolation between the PBE0 and the PBE. Many values of ω has been proposed in literature. Originally, the choice for the values of this quantity has been based on molecular tests and leads to a good reproduction of solids lattice constant, bulk moduli and band gap [51,52]. Two versions of HSE were developed, the HSE03 [48] with two different screening parameters for the HF and the PBE parts, ω HF = 0.15/ √ 2 bohr -1 and ω PBE = 0.15×2 1/3 bohr -1 . This specific version showed some problems. For example, the total energies were substantially lower than those found with the unscreened functional, this problem being strongly related to the different screening parameters used in HSE03. Therefore, a second more efficient version, HSE06 [54], with ω HF = ω PBE has been proposed. The HSE06 with a = 0.25 and ω = 0.11 bohr -1 gives a better account on the enthalpies of formation, ionization energy and ionization potentials, conserving a good accuracy for band gaps and lattice constants in solids [54].

Plane wave basis set

Once the exchange and correlation functional is fixed, the KS equation could be solved. The first step toward the numerical solution of KS Eq. [1.16] is to choose a basis set to expand the electron wave functions. In the case where the periodic boundary conditions (PBC) are used, often one uses a plane wave (PW) basis set, which is widely applied for the treatment of many condensed matter systems. Due to the periodicity of the potential, the KS eigenstates become Bloch functions and the single particle orbitals can be expanded as PW's. Using the Bloch theorem, the periodic wave functions of the systems could be written:

ψ k i ( - → r ) = e i( - → k . - → r ) φ k i ( - → r ) (1.28)
where ψ is the wave function of the periodic system, -→ k is the reciprocal space vector and φ is a function with the same periodicity as the system. Since φ is an arbitrary function, we could expand it as a PW's using Fourier series:

φ k i ( - → r ) = ∑ G C k i ( - → G )e i( - → G . - → r ) (1.29)
where -→ k is the Brillouin zone vector. The wave function reads:

ψ k i ( - → r ) = ∑ G C k i ( - → G )e i( - → G + - → k ). - → r , (1.30)
Where C k i is the Fourier coefficient and -→ G is a reciprocal space vector. The number of -→ G vectors can be tuned by making sure that the convergence of some basic properties is reached and considering only Fourier components corresponding to energies less than a cut off energy

E cut = 1 2 ( - → k + - → G ) 2 .
The choice of this cutoff energy depends on the studied system as well as the pseudopotential used (e.g next section) and determines the accuracy of the DFT energy.

Density Functional Theory

It is possible also to further simplify equation [1.30] by reducing the summation over k, for instance by including only the Γ (k=(0,0,0)) point. Such a choice is suitable for isolated systems and in general for cases in which the dispersion of the band structure can be safely approximated with a straight line. This is of course inaccurate for metals unless a very large unit cell is employed. Finally the wave function is given by:

ψ i ( - → r ) = ∑ G C i ( - → G )e i - → G . - → r . (1.31)
Plane waves are not explicitly dependent on the atomic coordinates. In addition, the accuracy of the basis set could be simply improved by increasing the number of the Fourier components used. However, as stated previously, the choice of E cut is always crucial. The number of the plane waves for a given E cut depends on the cell volume (Ω) and the K-points used. An estimation of the number of plane waves corresponding to a given cutoff is given by:

N PW = 1 2Π 2 ΩE 3/2 cut (1.32)

Pseudopotentials

Practical implementation of DFT encounter some problems. Any system is made of ions and electrons. The latter could be split into two kinds: the valence and the core electrons, the real contribution of these two type of electrons to the physically interesting properties of a given system being well established. It has been showed, since the early work of Fermi [18], that the action of the chemically inert core states exerted on the chemically active valence states could be replaced by an effective pseudopotential. This will lead to a reduction of the number of electrons to be explicitly considered and lowers the computational effort. An important advantage of pseudopotentials (PP) is the considerable reduction of the basis set size due to their smoother character compared to the core states. In addition, since the pseudized valence wave functions are nodeless functions, they allow to eliminate the rapid oscillations of the valence electron wave function in the core region. Figure [1.3] illustrates the above described image. In the core region, the pseudo wave function is nodeless and the pseudopotential does not have the asymptotic behavior of the all-electron one [56][57][58]. Outside the core region, the wave function and the pseudopotential match the allelectron ones.

Useful PP should take into account the core-electrons ion interaction, core-electrons Hartree potential energy, exchange correlation energy part and if needed the relativistic effects. From a mathematical point of view, the pseudopotential should be additive, i.e. the PP is the sum over all their respective pseudopotential contributions. Another crucial feature that a PP should satisfy is the transferability character, this property being required to allow the use of the same PP in different chemical environments. Among the large family of pseudopotential, the norm-conserving pseudopotentials is one of the most important. As it is called "norm-conserving", it consists on conserving the integrated density of the wave function and the pseudo wave function inside the radius r c [57,59,60].

Beyond the core region, the pseudo wave function should match as closely as possible the allelectron one.

∞ 0 |φ AE (r)| 2 r 2 dr = ∞ 0 |φ PS (r)| 2 r 2 dr (1.33)
Many authors proposed different types of norm-conserving PP such as Vanderbilt [61], Troullier and Martins [62,63], Kerker [64], Haman, Schluter and Chiang [59,65], Goedecker, Teter and Hutter [66]. The most general form of a pseudopotential reads:

V PP (r, r ′ ) = ∞ ∑ l=1 l ∑ m=-l Y * lm (ω)V l (r)δ (r -r ′ )Y lm (ω ′ ) (1.34)

Density Functional Theory

where l and m are the angular momentum quantum numbers and Y lm are the spherical harmonics.

The dependence of the PP on the coordinate r means that it has a local character. For each radial dependence, V l (r), we associate an angular momentum channel l described by its associated spherical harmonic Y lm (ω), where ω is the Euler angle of the position vector r which give it a non-local form. This representation is called semi-local form. The V l (r) can be reduced to the Coulomb potential (changing as -1/r for r → ∞) and as consequence it become fully radial. This allow to rewrite Eq. [1.34] as a sum of local and non-local PP contributions, where a suitable choice is made for L max :

V PP (r, r ′ ) = V loc PP (r) + L max ∑ l=1 l ∑ m=-l Y * lm (ω)∆V l (r)δ (r -r ′ )Y lm (ω ′ ) (1.35)

Total energy

The total energy derived from density functional theory in the framework of the plane waves basis set and the pseudopotential approximation could be derived from Eq.[1.13] by substituting the potential expression by the pseudopotential contribution of the local and non-local parts:

E KS [n(r)] = ∑ i ψ * i (r)(- 1 2 ▽ 2 +V ps nloc )ψ i (r)dr T KS +E H [n(r)] + E xc [n(r)] + E loc ion + ∑ I<J Z * I Z * J |R I -R J | E ext (1.36
) where

E loc ion = V ps loc n(r)dr (1.37)
The Z * I and Z * J are the total ions charge after removing the core electrons contributions. The associated Schrödinger like equation is given by:

(- ∇ 2 2 +V nloc +V e f f (r) H KS )ψ i (r) = ε i ψ i (r), (1.38)
The new effective potential reads:

V e f f (r) = V H +V loc ion +V xc (1.39)

Born-Oppenheimer molecular dynamics

DFT can be combined to the molecular dynamics methods to get a realistic description of the studied systems. This class of simulations go under the name of first-principles molecular dynamics (FPMD) simulations.

A priori, one can implement FPMD by using the so called Born-Oppenheimer molecular dynamics (BOMD). This method is based on the large difference in mass between the atomic nuclei and the electrons, the ratio being smaller than 10 -2 . One can assume that the electronic fast motion can be separated from the slow nuclei one, and that the electrons adapt instantaneously to the change on the nuclear configuration. Therefore, the electrons are taken as being always in their quantum mechanical ground state. BOMD consists on solving the electronic structure problem at each molecular dynamics step, given a fixed set of nuclear positions {R I }.

In other words, the electronic structure is first obtained by solving the time-independent stationary Schrödinger equation and then by propagating the nuclei's according to the classical dynamics laws in the potential obtained by solving the first step. The Born Oppenheimer molecular dynamics scheme is governed by the following equations:

M I RI (t) = -∇ I min ψ 0 {< ψ 0 [H e |ψ 0 >} (1.40) = -∇ I V BO e f f (R I ) (1.41
)

H e ψ 0 = E 0 ψ 0 (1.42)
where ψ 0 are the ground state wave functions and M I is the I th nuclei mass. Eq. [1.42] is the same as Eq. [1.38] derived for the ground state. As one can see from this set of coupled equations, at each step of the BO dynamics, the minimum of < H e > should be reached by diagonalizing the hamiltonian. Calculations of this kind are demanding, since at each step one has to bring the electronic system into its ground state, and perform a self-consistent field calculation. These difficulties may limit the simulated time scale and call for an alternative approach: the so-called Car-Parrinello molecular dynamics.

Equation of motion

The main idea of the Car-Parrinello method is to consider the electronic states {ψ i } as classical dynamical variables [67]. In this context the fictitious dynamical optimization of {ψ i } and the real atomic dynamics are run in parallel. When the nuclei move from one configuration to another, the electronic states {ψ i } are automatically optimized to the new configuration.

The fictitious dynamics of Car-Parrinello can be derived from the modified following Lagrangian:

L CP = 1 2 ∑ I M I Ṙ2 I + 1 2 µ ψ * i (r) ψi (r)dr -E({ψ i }, {R I }) + ∑ i j λ i j ( µψ * i (r)ψ j (r)dr -δ i j )
(1.43) The first term is the ions kinetic energy with M I the mass of the I th nucleus. The second term is the kinetic energy of the electronic variables with a fictitious mass µ. The third term is the total DFT energy. λ i, j in the last term is a set of Lagrangian multipliers introduced to keep the orthonormality of the Kohn-Sham orbitals. The Lagrangian determines the dynamic equations of the two set of degree of freedom:

M I RI = - ∂ E({ψ i }, {R I }) ∂ R I (1.44) µ ψi = - ∂ E({ψ i }, {R I }) ∂ ψ * i + ∑ j λ i j ψ j (r) (1.45)
Equation [1.44] describes the ions evolution on time at a temperature ∝ M I Ṙ2 I . The fictitious electronic motion is governed by Eq [1.45] and similarly evolve at a "fictitious temperature" ∝ µ ψ * i (r) ψi (r)dr. In this sense, if the fictitious electronic temperature is low enough, and if we optimize the initial configuration wave functions to its ground state, the electronic subsystem will remain close to its instantaneous minimum energy (minimum of the BO energy surface) during the dynamical evolution of the ions. The total conserved CP energy is given as:

E cons = 1 2 µ ψ * i (r) ψi (r)dr + 1 2 ∑ I M I Ṙ2 I + E({ψ i }, {R I }) (1.46)

Adiabaticity control

The key to achieve a proper handling of the CP method is to ensure the adiabaticity between the two families of degrees of freedom in such a way that the fast electronic subsystem follows the slow nuclear motion adiabatically. In the CPMD framework, the adiabaticity is monitored via the fictitious mass parameter µ. The control of the adiabaticity in practice could be checked by computing the vibrational density of state of the two subsystems. The two families of frequencies have to be well separated. Thus the energy transfer between the "hot nuclei" and the "cold electrons" cannot take place on the simulation relevant time scale. The value of µ is chosen so that the {ψ i } stay as close as possible to the exact BO energy surface minimum.

In practice, µ cannot be set to zero because in this case the dynamical optimization of {ψ i } will be prevented. Also, a very small value of µ will cause high frequencies to come to play, thereby reducing considerably the integration time step. However, large values are also not appropriate since they will increase the deviation from the correct ground state and affect the whole dynamics. As a consequence, the good value of µ has to be selected as a compromise between small deviations from the BO surface and sufficiently long time steps, compatible with affordable lengths for the time trajectories.

It has been shown [68] that the intrinsic dynamics of the orbitals can be described as a superposition of oscillations whose oscillation is given by:

w i j = 2(ε i -ε j ) µ 1/2 (1.47)
were ε i is the eigenvalue of the i th unoccupied and ε j the j th occupied level corresponding to

H KS e φ i = ε i φ i (1.48)
from Eq. [1.47] one can get that the lowest possible electronic frequency is estimated to

w min e ∝ E gap µ 1/2 (1.49)
Where the E gap is the electronic energy difference between the lowest unoccupied "LUMO" and the highest occupied "HOMO" orbital. This relation shows that this frequency increases like the square root of the E gap increases and decreases as the square root of the fictitious mass

w max e ∝ E cut µ 1/2 (1.50)
where E cut is the largest kinetic energy in an expansion of the wave function in terms of plane wave basis set as defined in Sec [1.2.3].

Decreasing µ shifts the electronic spectrum upwards on the frequency scale but it leads also to an increase on the maximum phonon frequency. Since the adiabatic separation of the electronic and ionic degrees of freedom is ensured when w min e is much higher than the highest ionic frequency w max n , the CP method should work well for systems with a defined gap. The value of µ should be fixed depending on the system. For a fixed µ the maximum integration time step is given by:

∆t max ∝ µ E cut 1/2 (1.51)
In case of metallic systems the use of the CP method is more troublesome, since they have a vanishing gap energy, so very small values of µ should be used to obtain the adiabatic separation. Indeed, the mere application of the Car-Parrinello ideas is bound to fail for metals, since, first and foremost, at finite temperatures, the number of electronic degrees of freedom depends on the temperature. For this specific kind of bonding, one has to refer to the free-energy molecular dynamics (FEMD) method [69] for a sound integration scheme accounting for finite electronic temperatures.

Forces in CPMD

The Car-Parrinello method consist on propagating the "hot ions" and the "cold electronic wave functions" simultaneously. For this purpose we need to evaluate the forces on each ion at each time step. The force should take into account the evolution of the Kohn-Sham eigen value when it moves to the new position. Thus, the forces acting on the orbitals are calculated as the action of the KS Hamiltonian:

F(φ i ) = -f i H KS |φ i > (1.52)
The forces acting on the nuclear are only the first KS energy derivative with respect to the ion positions R I :

F(R I ) = -∇ R I ψ i | Ĥ|ψ i = -(∇ R I ψ i |) Ĥ|ψ i -ψ i |∇ R I Ĥ|ψ i -ψ i | Ĥ(∇ R I |ψ i ) (1.53)
However, further simplification can be used if:

Ĥ|ψ i = E|ψ i (1.54) ψ i |ψ j = δ i j (1.55)
Hence, equation [1.53] is simply reduced to:

F(R I ) = -∇ R I ψ i | Ĥ|ψ i = -ψ i |∇ R I Ĥ|ψ i (1.56)
This formula is known as the Hellmann-Feynman theorem [70,[START_REF] Hellmann | Zur rolle der kinetischen elektronenenergie für die zwischenatomaren kräfte[END_REF]. Since the total energy is obtained from the DFT algorithms, the forces could be easily derived and thus the fictitious electronic and ionic motion are determined.

Thermodynamic control methods

Temperature control

Standard MD simulations are performed in the microcanonical ensemble (NVE) where the temperature is not a thermodynamic variable. Simulations at finite temperature in the NVT ensemble are also possible. To evaluate the temperature as an average over time, one needs to find a microcanonical observable for which the expectation value is a simple function of the temperature. In case of classical systems, this microscopic variable is the kinetic energy. Hence, the instantaneous kinetic temperature function for a system of N particle having masses M I and moving with velocities v i is simply given by:

N ∑ i=1 M i v 2 i 2 = 3 2 Nk B T (1.57)
The average temperature over the system over a trajectory made of M MD steps reads:

T = 1 M M ∑ m=1       1 3Nk B N ∑ i=1 M i v 2 i T (t m )       (1.58) 33 1.

Thermodynamic control methods

Where we assume the M MD steps taken into account correspond to uncorrelated configurations.

Velocity scaling

From Eq. [1.57], one deduces that if we can control in some way the velocities then we can keep under control the related parameter T. If at a time t, the temperature is T (t) and the velocities are multiplied by a factor λ , then the temperature change with respect to the target temperature T 0 can be calculated as

∆T = T (t) -T 0 = 1 3Nk B N ∑ i=1 M i (λ v i ) 2 - 1 3Nk B N ∑ i=1 M i v 2 i (1.59) ∆T = (λ 2 -1)T 0 (1.60) λ = T (t)/T 0 (1.61)
The temperature control consists in multiplying the velocity at each time step by λ . This method allows to control the instantaneous temperature, but it does not correspond to the practical realization of a the canonical ensemble (see below).

Nosé-Hoover thermostats

A rigorous method to control the temperature was originally proposed by Nosé [72,73] and subsequently developed by Hoover [74]. The original idea of Nosé consists in adding a new artificial variable s, associated with a velocity ṡ and a fictitious "mass" Q > 0. The parameter Q has the dimension of energy time square and determines the coupling between the heat reservoir and the real system and so the time scale of the temperature fluctuation. The artificial new variable s is a time-scaling parameter that ensures the canonical sampling. Similarly, for the electronic degrees of freedom; the thermostat extension was proposed by Blöchl and Parrinello [75], by adding a single new variable ξ associated with a fictitious mass Q e . The rescaling of the electronic and ionic velocities can be obtained as follows:

ṘI → s. ṘI (1.62) ψi (x) → ξ . ψi (x) (1.63)
The virtual potential for s(t) is given by the Boltzmann-like canonical term (3N + 1)k B T lns, where 3N + 1 are all the degrees of freedom: 3N for the atomic coordinates plus one for the new variable s. The virtual potential for ξ (t) is given as -2E kin e lnξ . E kin e is the CP fictitious electronic kinetic energy at which we want to keep the electrons.

The choice of this parameter will be discussed later. Then, the extended CP Lagrangian with the new variables reads:

L =µ ∑ i ξ 2 | ψi (x)| 2 + 1 2 ∑ I M I s 2 Ṙ2 I -E DFT (1.64) + ∑ i j λ i j ( ψ i |ψ j -δ i, j ) (1.65)
+ 1 2 Q e ξ 2 -2E kin e lnξ (1.66) + 1 2 Q ṡ2 -(3N + 1)k B T lns (1.67)
The coupled Euler-Lagrange equations of motion, under the rescaling of the Nosé-Hoover thermostat become then:

µ ψi (x) = - δ E DFT δ ψ * i + ∑ j λ i j ψ j (x) - ξ ξ .µ.ψ i (x) (1.68) M I RI = 1 s 2 ∇ R I E DFT -2 ṡ s M I ṘI (1.69) Q e ξ = 2 ξ ξ 2 ∑ i µ| ψi (x)| 2 -E kin e (1.70) Q s = 1 s s 2 ∑ I M I Ṙ2 I -(3N + 1)k B T (1.71)
If we check equation [1.69], the term 2 ṡ s can be considered as a dynamical friction coefficient. This leads to a dissipative dynamics with a non-Hamiltonian behavior. The sign of the friction term could be either positive or negative depending on its equation of motion [1.71]. This will generate a kind of acceleration or damping of the ions and thus a heating or cooling of the subsystem. As a consequence and in contrast to the case of the velocity scaling the extended dynamics based on Nosé-Hoover thermostat control can be shown to reproduce the canonical 34 1.5. Thermodynamic control methods ensemble fluctuations in the nuclear and momentum subspace. Similarly, the electronic dynamics show the same canonical behaviour. A guideline value of the fictitious kinetic energy of the electronic subsystem could be estimated under the assumption of a canonical model of a well separated atoms [75]. In this model the wave functions follow the corresponding atom rigidly so we can relate the fictitious kinetic energy to the atom velocities. If we assume that all the atoms have the same mass, the fictitious kinetic energy reads:

E kin = 2k B T µ M ∑ i ψ i | - 1 2 ∇ 2 |ψ i (1.72)
Where M is the mass of one atom. The optimum value recommended for practical applications is about twice E kin given by eq. [1.72].

Pressure control

The necessity of a theoretical counterpart to the constant pressure and temperature experiments require the development of a theoretical framework compatible with first-principles molecular dynamics simulations. The most natural simulation technique to mimick the application of pressure consists in working in the NVT ensemble with a larger reduced density. This method demands the knowledge of the ambient pressure system density ρ 0 , then the new system under pressure will be generated at a new density ρ > ρ 0 by changing the atomic coordinates as a result of a reduction of the volume. Despite the simplicity of this method, the target pressure can not be estimated beforehand unless the equation of state is well known. The pressure can be computed either during the simulation or at the end of the simulation after relaxing the system structure to zero temperature.

In principle, one can also calculate the pressure (third of the trace of the stress tensor) and calibrate the volume by trial and error. This method is cumbersome and can require several adjustments. Another way to control the target pressure is to employ the NPT ensemble. The application of an external constant pressure leads to volume fluctuations. Thus, the volume should be considered as a new dynamical variable. A number of different barostat techniques exist for maintaining a target pressure by way of adjusting the simulation volume. In the following we review the Anderson barostat since this is the scheme employed within this work.

Andersen barostat

In the method proposed by Anderson [76], the volume is allowed to change only for an isotropic variation where the cell shape remains unchanged. In this case, the volume Ω = l × l × l is considered as a new Lagrangean variable were l is the box side length:

L = 1 2 N ∑ I=1 M I l 2 ṙ2 I -V (lr I ) + 1 2 Q Ω2 I -P 0 Ω (1.73)
Where P 0 is the external pressure, Q is the fictitious mass parameter associated to the volume dynamics. For convenience, the atomic coordinates are rewritten in a scaled form as:

r I = R I /l (1.74)
the new r I variable is not related to the volume which is henceforth a completely independent dynamical variable. From Eq.[1.73] one could derive the equations of motion for the particles r I and the volume Ω which read:

M I rI = 1 l f I - 2 3 Ω Ω M I ṙI (1.75) Q Ω = P -P 0 (1.76) 1.6. CPMD implementation
where f I is the force from Eq. [1.1]. The quantity P is the dynamical computed pressure given by:

P = 1 3Ω ∑ I M I l 2 ṙ2 I + ∑ I>J (R I -R J ) f IJ (1.77)
The only new input parameters are the fictitious mass parameter Q and the target external pressure P 0 .

CPMD implementation 1.6.1 Wave function optimization

The wave function optimization is a prerequisite to any CP molecular dynamics. It consists on getting the electronic system into its ground state for the initial ionic configuration. In the BO molecular dynamics scheme, this step is embedded in the dynamical equation, while in case of CPMD the Schrodinger equation is not solved at each step, this is why we need to make sure that at the beginning the electronic ground state has been properly obtained. The CPMD keywords needed are displayed in this table: &CPMD . . .

OPTIMIZE WAVEFUNCTIONS INITIALIZE WAVEFUNCTIONS RANDOM TIMESTEP

5 MAXSTEP 1000 EMASS 1000 PCG CONVERGENCE ORBITALS 1.0E-06 . . .

&END

The wave function could be either initialized from a random choice or from atomic positions, then the keyword AT OMS is used instead of RANDOM TIMESTEP: is the time step used to integrate the equations of motions MAXSTEP: is the maximum number of the optimization steps, default value is 10000 EMASS: is the electronic fictitious mass parameter PCG: is the preconditioned conjugate gradients method for the optimization CONVERGENCE ORBITALS: is the convergence criteria, the default value is 10 -5

A typical CPMD output is the following: The columns have the following meaning: NFI:

NFI GEMAX CNORM ETOT
Step number (number of finite iterations) GEMAX: largest derivative of the off-diagonal Fourier component CNORM: average square of the off-diagonal Fourier components ETOT: total energy DETOT: change in total energy to the previous step TCPU: (CPU) time for this step The ground state is reached one the GEMAX is lower than the convergence criteria, in the above example: 10 -6 a.u.

Car-Parrinello Vs Born-Oppenheimer conserved quantities

In order to compare the CPMD and the BOMD schemes, both methods were tested with different time step parameters and different convergence criteria in the case of a GeSe 2 glass at 38 in particular the system total energy as computed from Eq. [1.46].

The top panel features the comparison between CP1 and CP2. Both dynamics show a very good stability, the best energy accuracy being associated to the CPMD run with an time step of 5 a.u.

In this case we have an energy conservation of 5 × 10 -9 a.u/ps. A larger time step (10 a.u) (eg. CP2 run) leads to a larger energy fluctuations with an energy conservation of 1 × 10 -8 a.u/ps. The middle panel displays the CP2 run (black line), BO1 run (green line) and BO3 (magenta line). The CP run gives a more precise energy conservation than any of the BO run. For this particular system, the BO1 run exhibits a large inaccuracy compared to the stable CP2 result. Since the energy conservation decreases drastically to 2.6 × 10 -4 a.u./ps. Surprisingly, when increasing the time step to 100 a.u., we recover the CP accuracy with an energy conservation of about 6 × 10 -6 a.u./ps. This behavior remains to be explained. Since decreasing the convergence criteria leads to a lower computational effort, we tried a run within the BO scheme and with a convergence of 10

CPMD implementation

Overall, the CPMD schemes leads to a very good energy conservation, while the most stable BOMD run was performed with a time step of 100 a.u. and a convergence of 10 -6 . One may ask if the CPMD performances remain unchanged when high temperature are reached. The following figure show an example of simulation of an amorphous GeSe 2 for a complete thermal cycle. Therefore, the high temperature results are also very stable.

Turning now to the computational performances, the time needed to generate one picoseconds of Car-Parrinello MD trajectory decrease with increasing the time step, similarly to the Born-Oppenheimer case when the convergence criteria is kept the same. BOMD speed up as expected when decreasing the convergence, with a worse energy conservation. For a general instructive view on the CPMD code, we performed tests on different platforms and with different system sizes. Figure. [1.7] The overall remark is that CPMD goes faster when increasing the number of processors used even with large number of CPUs (fig (c)). The use of many OMP threads may slow the calculation speed because it increase the communication time between the cores. CPMD show a good scalability on different platforms and a high level of parallelism.

Velocity scaling

The equation that define how to control temperature via velocity scaling are detailed in Sec.[1.5.

&END

For instance, we imposed to the ions to thermalize at T=900 K with a window of plus or minus T=50 K. For the electronic degrees of freedom, the target kinetic energy is of 0.04 a.u. with an acceptable variation over the range of plus and minus 0.01 (a.u.). An illustrative example is displayed in the following picture: This method is not recommended since it amounts to subtracting energy from the system so as to create discontinuities in the energy variations. Also the canonical sampling can not be reproduced. However, for a some particular cases this method could be a practical solution (e.g high temperature).

Thermostat frequencies

As explained in detail in section [1.5.1.2], the temperature control is ensured by the use of the Nosé-Hoover thermostats that control both ionic and electronic degrees of freedom temperatures. From equations [1.70] and [1.71] it appears that in addition to the ionic temperature T and the electronic kinetic energy E kin , the fictitious masses Q and Q e are also input parameters.

The choice of these parameters will impact the accuracy of the thermostat control.

Once the target ionic temperature is fixed one could estimate the adiabatic fictitious kinetic energy E kin from the Blöchl formula Eq. [1.72]. In practice, the CPMD keywords needed to keep a good thermodynamical control of a GeSe 2 system are presented in the following table:

&END

Here we fix the ionic temperature (T=300 K) and the E kin (0.04 a.u.) and, instead of the fictitious masses associated to each thermostat, we give the correspondent frequencies. The frequencies, ω R associated to the ionic subsystem thermal control can be estimated in the framework of the harmonic oscillator relationship as:

ω R = 2gk B T Q (1.78)
Similarly, the frequencies associated to the thermal fluctuation of the thermostat that control the electronic degree of freedom, can be estimated as follows:

ω e = 4E kin Q e (1.79) 1.6. CPMD implementation
As the choice of the fictitious electronic mass is of importance to achieve the adiabatic conditions, the choice of the thermostat frequencies has the same effect. Thus, in order to get an adiabatic well separated dynamics of the different families of degree of freedom, some rules should be respected. Figure [1.9] illustrate the frequencies scale. In order to recover the adiabatic separation, the electronic frequencies have to be higher than the ionic frequencies in such a way there is no any overlap between the two sets of frequencies. The frequencies of the "electronic" subsystem thermostat have to be lower than the frequencies of the "electronic" subsystem. Similarly, the frequencies of the ionic thermostat subsystem should also be kept lower than the frequencies of the ionic subsystem. The numeric values of the ionic thermostat frequencies have to be fixed at first, for example by computing the dimer vibrational frequency, then the electronic thermostat frequencies is "chosen" to ensure a large frequencies gap between ions and electronic degrees of freedom. For example, let's consider a system of GeSe 2 with 120 atoms (40 Ge and 80 Se atoms) at T= 900K, by applying the Blöchl formula Eq. [1.72]. The kinetic energy was found to be E kin = 0.012 (a.u). Thus, as recommended, we fixed the target kinetic energy to E kin = 0.04 (a.u). Figure [1.10] shows ionic and the electronic kinetic energy during a typical MD run at different temperature. The ionic thermostat frequency was fixed to 200 (a.u) and the electronic fictitious kinetic energy to 0.04 (a.u) with a frequency of 600 (a.u). The above figure shows that the thermostat ensures a good control of the kinetic energy (or the temperature) for all temperatures and during a long simulation period. Unfortunately this good control is not always manageable if the system is subject to extreme thermodynamical conditions. Figure [1.10] shows an MD run on the same GeSe 2 system with the same thermostats parameters but with an additional constrain. In fact the system is subject to a pressure of 9.8 GPa.

At high temperature, while the ionic thermostat ensure a good control of the temperature, the electronic thermostat does not perform as well. In particular, a departure from a constant behavior is observed between t=64 ps to t=115 ps. This behavior can be ascribed to a change in the bonding nature of the system or, more simply, to some level-crossing of the electronic states, bringing into the dynamical evolution new variable for which no previous dynamical evolution existed.

Barostat fictitious mass

Maintaining a constant pressure using the Parrinello-Rahman method (or its version limited to changes of volume size, with a fixed shape) requires a proper definition of the fictitious mass parameter ensuring a good stability of the cell dynamics (see Eq. [6.60], appendix E). By referring to Ref. [77], a value of W were optimized for a NiZr 2 crystalline system simulation in the NPT ensemble and it was set to W * = 8.61 × 10 7 (a.u). I took this value as our initial one to optimize the barostat fictitious mass W for the systems of our interest. CPMD typical keywords for an NPT simulation are illustrated in the following figure:

&END

For illustration, many values were tested for amorphous GeSe 2 at T=900 . For this test, a high temperature is recommended because at ambient and low temperature all the tested values gave a good stability of the conserved quantity E tot , while at high temperature some difficulties in the choice of the optimal value may appear. Figure [1.11] display the results of our tests. For W = 10 6 a.u. the dynamics is highly dissipative and the conserved quantity is not anymore conserved. This problem is less severe for W = 10 7 a.u whereas for values higher than 10 8 (a.u) an equilibrium value. For higher values (W = 10 7 a.u. and W = 10 8 a.u.) the cell variation is smoother. However, for very high values (W = 10 9 a.u. and W = 10 10 a.u.) the mass parameter becomes very heavy and leads to very slow cell dynamics. Roughly speaking, this is similar to a cell not following the dynamics of the system, that is to say, a volume that remains essentially constant as in the microcanonical ensemble. Based on this tests, the value that guarantees a good stability of the total energy and, at the same time, a good cell dynamics is chosen to be W = 10 8 a.u.. Henceforth, all our NPT simulations have been carried out with this value.

a

Calculating the Pressure

The pressure is computed from the total internal stress tensor formula Eq.[ In this case the total stress is equal to 0.27 GPa. One should note that for all our NVT calculation the stress tensor can been computed at T=0 K after annealing the structure and taking out all the thermal contributions. For the NPT simulations, the values that we display were computed by taking the internal stress average over the configurations collected at T=300 K. The following figure displays the average values of the diagonal term (equal to the internal stress) and the off diagonal terms of the stress tensor for an amorphous GeSe 2 simulation in the NPT ensemble at ambient conditions: The average computed pressure over the 100ps of dynamics is 0.13GPa, different from the target value P 0 = 0GPa. In order to understand this difference, we have calculated the kinetic part of the pressure given as an average by:

P T hermal = NK B T V (1.80)
Numerical application to the system displayed in the figure [1.13] with N=120, V=3679.0289 Å 3 , K B = 1.380610 -23 m 2 Kgs -2 K -1 and T= 300K, knowing that 1GPa = 10 9 Kgs -2 m -1 one find a thermal contribution to the stress P T hermal = 0.13GPa. This results points toward a possible inconsistency in the coding of Eq. [1.76], since the value that we calculated is essentially the same as the kinetic part. All this is consistent with an erroneous expression of P in Eq. [1.76] based on the potential part only. This part is set to zero (in average) with no consideration of the thermal part. For these reasons P = P T hermal = 0.13GPa.

the targeted systems, our focus will be on to the Ge x Se 1-x family, in which the nature of the network arrangement and the bonding strongly depends on the composition x. Within this family, glassy materials form easily for x < 0.43, while for x = 0 the network is mainly made of Se rings and chains. As a guideline, phenomenological ideas based on rigidity can be exploited to achieve a first picture of these networks. Indeed, the rigidity of the network is a well established Taken from [1].

property, that depends on the number of constraints due to the bonds in the network. Based on this criterion, three concentration regions with a special network connectivity can be identified Fig. [2.1]. The under-constrained network for x < 0.2, the over-constrained network for x > 0.25, while in between the two regions an intermediate window exists where the network is optimally constrained in the sense that the number of constraints is equal to the number of degrees of freedom of the network. This transition region is marked by a floppy network. Over these ranges, the physico-chemical properties of the glass and liquid materials exhibit drastic structural and electronic changes. Furthermore, the physical origin and the manifestation of these changes is not always obvious to track under extreme constraints. However some similarities and marked signs do exist in all materials. In the following, we will present the main results and historical steps toward the actual understanding of the atomic scale properties of binary chalcogenides.

GeSe 2

The first "DFT-inspired" molecular dynamics results on binary GeSe 2 were obtained in the middle ninetieth by Cappelletti and Co. [2]. The authors studied at first the vibrational properties of the glassy state and compared them to neutron scattering experiments. The scheme used is the one developed by Sankey and Co. [3] based on a multicenter tight-binding model for molecular dynamics founded on density-functional theory. The tight binding model is a method used to compute the electronic band structure using an approximate set of wave functions that are assumed to be a superposition of wave functions of isolated atoms located at each atomic site. The DFT approximated energy is based on the Harris [4] and Foulkes [5] approximation, for which a reference electronic density is taken as a sum of neutral-atom spherical atomic density. By using this approximation, an energy functional could be expressed based, for instance, on the LDA exchange-correlation recipe. It should be reminded that this scheme is a non self-consistent one, since the electronic density of the system is a superposition of isolated electronic densities of atoms or weakly coupled fragments.

The vibrational properties calculated via the Harris functional approximation are in good agreement with the measured one. Despite the fact that the model was very small (63 atoms), it reproduced successfully the regions of the A 1 and A 1c bands characteristic of the GeSe 2 glass whereas a significant difference occurred near 20 meV. The computed structure factor is in good agreement with the neutron-scattering results (see Fig. points are experimental results from ref [6]. Figure taken from [2].

Overall, this work does not give any explanation either of the origin of the vibrational modes either of the first-sharp diffraction peak (FSDP) observed in the computed static structure factor. A second similar study [7,8] within the same theoretical framework and with a larger model (216 atoms), gave improved structural and vibrational properties, also providing further structural analyses. The need of an accurate quantitative description of the structure of binary alloys require to go a step further toward realistic simulations based on a fully self consistent theoretical framework. The first set of fully self consistent first-principles molecular dynamics calculations on a disordered GeSe 2 system was published by Massobrio and Co. [9,10] on the structure of the liquid GeSe 2 . The first issue addressed was the quantitative character of the model employed for chalcogenides. In particular, these authors expressed some concern on the choice of the exchange and correlation functional to be used. Although at that time this set of results were not immediately published, the first calculations in this respect were carried out using the local density approximation (LDA) for the exchange and correlation part of the DFT total energy. These calculations revealed that the LDA results were highly disappointing and unable to reproduce essential features of the experimental data, such as the FSDP in the total neutron structure factor. Years later and for the sake of consistency, the LDA results were shown and compared to those obtained by using the generalized gradient approximations for the exchange-correlations functional (GGA, in particular the Perdew-Wang (PW) exchange-correlation one [11]). As shown in figure [2.3] the LDA GeSe 2 liquid does not feature the first sharp diffraction peak (FSDP) in any of the Faber-Ziman (FZ) structure factors.

In order to understand why such very different structural results occur when two different XC functional are used, as for instance the Perdew Wang (as a case of GGA) and the LDA functional, a detailed rationale was put forth by Massobrio and Co. [13]. In this work the effect of changes in the electronic structure and how this affects the chemical bonding (as an electronic property) and the intermediate range order (as a structural property) were highlighted. Since the LDA and the PW-GGA approximations account differently for the exchange and correlation energy, different degrees of ionicity (and, in principle, of covalency) had been automatically introduced by performing distinct calculations in the two separate frameworks. Also, information is gathered by optimizing the electronic structure (separately via both the LDA and PW-GGA schemes) for specific configurations obtained within one given framework. To correlate the changes in the bonding character to a given functional, the electronic structure and the valence charge distribution had been analyzed for selected configuration of both models by forming a specific Se-Ge-Se trimer (see Fig. Ref. [13] The investigation of the valence charge profiles revealed an accumulation of charge around the Se atom and a depletion around the Ge atom resulting from the higher electronegativity of Selenium. In addition, the existence of lobes pointing along the bond directions indicates the importance of the covalent contribution to the bonding character. When deeply looking at the counter plots, the GGA seems to favor a charge transfer compared to LDA, thereby pointing out the higher degree of ionicity. This is confirmed by the contour plot of the difference between the charge densities obtained with the two functionals. The spherical nature of the charge distribution in the difference contour plot has important implications. The distinction between the two schemes lies in the way the ionic character is described. The PW-GGA favor a higher degree of ionicity, the covalency remaining equivalent in the two schemes. As a consequence, the occurrence of the FSDP in the neutron structure factor can unambiguously be related to an increase in the ionic character of the bonds, since the LDA drastically fails to reproduce its intensity and position, Fig. [14]. For clarity, LDA curve is displaced downward by 0.4. From ref. [15].

In the quoted paper, it was emphasised that the LDA scheme compare only qualitatively with experiment, being inadequate to describe the short-and intermediate-range order. For this reasons the PW-GGA was taken as the method of choice to generate, study and unravel the structural properties of binary chalcogenides. The first applications of this theoretical recipe was the investigation of the structural properties of liquid GeSe 2 [9]. In these work, the authors gave an atomic scale description of the GeSe 2 liquid state of unprecedented quality. The first focus was on the reproducibility of the FSDP characterizing the total structure factor measured experimentally [12] and calculated in classical MD [16]. The second issue was related to the occurrence of the FSDP in the Bhatia-Thornton concentration-concentration structure factor (S cc (k)). Experimentally [12], this peak was clearly present, as opposed to what found in the study based on empirical potentials [16][17][18] where it was conjectured that the absence of the FSDP in the S cc (k) structure factor could be a general property generic to all binary glasses. The third issue addressed was related to the chemical order. Indeed, the classical MD studies revealed that the GeSe 2 is mainly an ordered network of interlinked tetrahedra. This is contrasted by experimental evidence, showing that a departure from the chemical order occurs in the liquid and the amorphous states [12]. FPMD-DFT calculations were performed in the NVT ensemble. In this approach, the valence electrons are treated explicitly, whereas norm-conserving pseudopotentials are used to account for the core-valence interactions. As stated above, the exchange and correlation functional due to Perdew and Wang (PW) was used. The simulated system contains 120 atoms in a periodically repeated cubic cell of size 15.7 Å, corresponding to the liquid density at 1050 K. The GeSe 2 network was found to be more subtle than the theoretical descriptions provided by the continuous random network (CRN) and the crystalline like models. The Ge atoms do exist in different environments, 63% of Ge atoms in a four fold environments coexist with 30% of non bonded ones. The Se atoms are found in two and three fold environments with respectively 74% and 23%. In contradiction with the classical MD simulation, a non negligible fraction of Ge (10%) and Se (39%) was found to exist in homopolar configurations confirming the departure from the chemical order. Turning to the reciprocal space properties, the total structure factor was nicely reproduced having a very good agreement with the measured one including the low wave vectors region where the FSDP is expected to occur Fig. experimental data (circles). From ref. [9].

When comparing the computed Faber-Ziman partial structure factors to the measured ones, the general shape is reproduced with a high accuracy, whereas the intensity of the FSDP is underestimated in the case of Ge -Ge correlations, and overestimated in the case of the Ge -Se correlations. In addition, the S cc (k) Bhatia-Thornton structure factor FSDP was still absent in line with Vashista's statement. As a corollary motivation for this work, emphasis was given to the origin of the intermediate range order (IRO) through the appearance of the first sharp diffraction peak (FSDP) [19].

The FSDP has a peculiar behavior with temperature, especially when the temperature is raised highly beyond the melting point [12,20]. To check this statement, attention was devoted to the correlation between the number of predominant tetrahedral units, the intensity of the FSDP and the temperature. By using the same computational scheme, a high temperature model was generated at 1373 K and was compared to the experimental counterpart one [12,20]. The comparison of S(k) with experiments shows a good agreement in the reproduction of the FSDP. In particular its intensity exhibits a huge reduction at high temperature (see Fig. Open circles: experimental results, from [12] (T = 1050 K) and [20] (T = 1373 K). Solid line: theory, [9] (T = 1050 K) and [19] (T = 1373 K). From Ref. [19] Further investigation showed that there are no cancellation effects since the partial structure factors exhibit very low intensities of their different FSDPs. The decrease of the FSDP was correlated to the reduction in the number of GeSe 4 tetrahedral conformations. Indeed, the fraction of four fold Ge centered motifs decreases, while three-fold Ge and miscoordinated Ge atoms appear instead. Also, the fraction of Se -Se homopolar links increases along with the temperature generating an important change in the Se local environment. This study showed that the breakdown of the IRO is due to the appearance of various chemically disordered motifs, leading to a distribution of structural units devoid of any predominant (higher than, say 40-50 %) motif.

As an additional contributions toward the understanding of the origins of the FSDP, it appeared of interest to correlate its intensity with the range of correlations accounted for in real space. To this end, the S FZ Ge-Ge structure factor has been calculated for an extended cell of 3240 atoms by Fourier transforming (FT) the real space pair correlation functions (PCF) [11]. This extended cell has been constructed by replicating the 120 atoms cell in the three directions. First, the maximum possible cut-off distance in real space was used to Fourier transform S FZ Ge-Ge , the comparison showing a pattern very close to the directly computed one. Second, by reducing the cutoff radius (r c < r max ) (e.g see Fig. [2.8]), it appeared that for r c larger than 10.5 Å a well defined peak still persists at the FSDP peak location. However, for values smaller than 6 Å only a residual bump remains at the FSDP position, while for 5 Å this feature disappears completely. As a consequence, it was concluded that the FSDP in the S GeGe (k) structure factor is mainly due to correlations beyond 5 Å . A similar analysis has been employed for the total structure factor, confirming that the range (r c > 5 Å ) is responsible of the appearance of the FSDP. For r c = 5 Å a shoulder persists between 1 Å -1 and 1.5 Å -1 which proves that the FSDP carries also some contributions from Ge -Se and Se -Se correlations.

Figure 2.8 S tot

neutron (k) of l -GeSe 2 (solid lines) obtained by Fourier integration of the calculated PDF within given integration ranges 0-r c . The dotted line corresponds to the S tot neutron (k) directly calculated in k space. From ref. [11] 2.2. GeSe 2 Many alternatives hypothesis were proposed to correlate the FSDP to specific structural motifs. Two famous suggestions were largely discussed, the first links the FSDP to the presence of a crystalline like signature where its position is related to the inter-layer separations [21,22]. The second claims the occurrence of cluster regions separated by interstitial voids giving rise to correlation distances typical of the intermediate range order [23]. The reliability of these hypothesis was tested by studying the liquid structures of the GeSe 2 and SiO 2 (the theoretical framework being unchanged, namely FPMD simulation [24]). This work brings direct evidences based on the calculation of the second order moment of the structure factor. It was found that the two frequently invoked interpretation schemes, based on the occurrence of layered structures and of cluster-void correlations, are not totally convincing to account for the appearance of a FSDP.

Overall, PW-GGA calculations on liquid GeSe 2 were instrumental in setting the scenes for a quantitative description of disordered chalcogenides based on FPMD. However, the quest for a real predictive power called for the search of improvements, very much needed especially when focusing on Ge-Ge correlations. To do this, the choice was to explore the effect of different exchange-correlation (XC) functionals. The GGA family of DFT exchange and correlation functional do contain different proposals, as stressed in the theoretical chapter [see section 1.2.2.2]. For instance, the generalized gradient approximation after Becke (B) for the exchange energy and Lee, Yang and Parr (LYP) for the correlation energy is an alternative to go beyond the PW description. The reliability of BLYP has been assessed in different DFT calculations. It has been shown that the BLYP functional performs very well for properties encompassing equilibrium geometries, vibrational frequencies, and atomization energies of nanostructured systems. The success of the BLYP is in part due to the LYP correlation energy, in particular to the correlation energy formula due to Colle and Salvetti [25] which was recast in terms of the electron density and of a suitable Hartree-Fock density matrix, thus providing a correlation energy and a correlation potential. This scheme is expected to enhance the localized behavior of the electron density at the expense of the electronic delocalization effects that favor the metallic character. In the case of liquid GeSe 2 , it has been shown that choosing BLYP improves upon PW and LDA results [26]. In particular BLYP minimizes the electronic delocalization effects that are unsuitable in any bonding situation characterized by competing ionic and covalent contributions. The simulations based on BLYP exchange-correlation functional have shown a substantial improvement of the short-range properties of liquid GeSe 2 . The shapes of the Ge -Ge and to a lesser extent the Ge -Se pair-correlation functional was more structured which is consistent with a higher level of tetrahedral organization (see Fig [2.9]). [11] , open circles: experimental results of Ref. [12]. From Ref. [26] The BLYP mean coordination number is only 3.5% higher than the experimental one instead of 7.7% in the PW case. The bond angle distributions (BADs) give a sharper distributions, the Se -Ge -Se is more symmetric around the tetrahedral angle 109 • while the Ge -Se -Ge shows two distinct peaks resulting from edge-sharing (ES) and corner sharing (CS) connections. The impact of the BLYP scheme on the intermediate-range properties remains elusive. Smaller improvements are found for the intensity and the position of the FSDP in the partial structure factor [27]. The impact of the BLYP XC functional has been also assessed through the modeling of the amorphous GeSe 2 [27,28]. A clear improvement has been found on the real space properties as well as on the reciprocal space properties with respect to the PW scheme. The Ge-Ge pair correlation function was found to be in a good agreement with the measured 2.3. GeSe 4 one. This agreement is substantiated by the reproduction on the homopolar bonds fractions with respect to the experimental one. In summary, the modeling of the prototypical network forming GeSe 2 system, briefly summarized in this section, can be considered as resulting from three different research efforts. First, interest has been given to choose the better XC functional and it was emphasised that the GGA family is suitable to improve the predictive power of the models. In addition, a good account of the ionicity through the use of GGA functionals was found to enhance the theoretical description. Second, interest has been devoted to practical application of the established theoretical scheme to study disordered chalcogenides. In particular, pieces of evidence have been provided for the high temperature effect, the origin and the reproducibility of the FSDP in the partials structure factors and in the total structure factor. Third and last, the focus was turned to the XC functionals within the GGA family where the BLYP proved to give better descriptions of fundamental structural properties, and reproduces more accurately delicate features of the liquid and glassy states.

GeSe 4

After the success of the FPMD methods to elucidate the microscopic structure of GeSe 2 interest has focused on GeSe 4 . In the work by Haye and Co. [10], the total structure factor showed a good agreement with the diffraction data for a similar concentration (x = 0.15). The existence of the FSDP indicates that the IRO in the liquid state is well described. The real space properties as described by the partial radial correlation function revealed a good agreement between the calculated partial coordination numbers and the 8 -N phenomenological rule. The obtained coordination numbers are also very close to those predicted by the Chemically Ordered Network model (CON). The Ge atoms are found to be mostly on a four fold conformations (86%) with a small fractions of three and five fold atoms. The Se atoms are found in a two fold geometries (89%) as predicted by the 8 -N rule and the CON model. The microscopic structure is described in terms of GeSe 4 tetrahedra that coexist with a high fraction of homopolar Se chains. Two extreme situations can be considered to describe the network. In the first, a phase separation exists between GeSe 2 and chains of Se atoms. In this case there are no Se atoms connected to one Se atoms and one Ge atom. Alternatively, one can think of Se atoms connected to Ge and Se atoms, thereby ruling out Se-Se-Se connections and Ge-Se-Ge connections. In this case the tetrahedra are fully interlinked by Se atoms. It appears that Se atoms can be classified in terms of the bonding types, namely Se-Se-Se, Se-Se-Ge or Ge-Se-Ge (e.g see Fig [2.10]). to the extremely high quench rate used by that team. The model of glassy GeSe 4 at the origin of paper [29] was used by Kibalchenko and Co. [31] to compute the density-functional NMR chemical shifts. The analysis compared the theoretical NMR spectrum to the measured one and derived a composition of the studied glass in terms of local coordination. This study is consistent with a decomposition in terms of structural motifs which coherently accounts for both NMR and neutron diffraction experiments. About 20% of the Se atoms are found in Se-Se-Ge linkages (AB type). For GeSe 4 in particular, the results were not consistent with the phase-separation model, supporting the occurrence of a fully bonded network in which Se-Se-Ge linkages ensure the connection between Se-Se-Se chains and Ge-Se-Ge tetrahedral arrangements.

The atomic structure of glassy GeSe 4 has been further studied in the framework of firstprinciples molecular dynamics [32], by considering another XC functional, namely the PBE (Perdew, Burke, Ernzerhof) one. It has been found that the short range structure in glassy GeSe 4 is only moderately modified when switching from the Perdew-Wang to the Perdew, Burke, Ernzerhof GGA scheme. For instance, small shifts in the positions for the g GeGe (r) pair correlation function were found. The impact of the PBE recipe is reflected in the g GeGe (r) pair correlation function, where we found a more intense peak and a deeper minimum. This is associated with the increase of the fraction of edge-sharing motifs (PBE:30%; PW:22%) (e.g see Fig. [2.11]).

Interestingly the PBE calculations highlight the presence of extended Se n chains (n = 12) which are compatible with a network structure in which GeSe 4 and Se n units of various lengths coexist in a fully interconnected fashion. Again, no sign of phase separation was found, this closing the debate on the possible occurrence of this phenomenon, found not to take place on the theoretical point of view for three different exchange-correlation functionals within DFT.

Motivation

As I sketched and summarized in the previous two sections, building realistic models and reliable methodology to study amorphous materials required a slowly step by step evolution.

Actually and after solving the main issues related to the FSDP and the IRO, interest should automatically be given to the possibility of exploring their phase diagram. In particular, the structural and the electronic properties changes that could happen when a given structure is subject to high pressure or temperature are worth of interest. In the first part of this thesis our interest is mainly devoted to high pressure effects. Experimentally, the crystallization and the amorphization of binary amorphous systems under pressure was studied. In 1965, Prewitt and Young [33] quoted that under pressure and at T=300 K two amorphous phases were found. Later, Shimada and Dachille [34] studied the crystallization of amorphous GeSe 2 and GeS 2 under pressure. They showed the existence of three GeSe 2 crystalline phases depending on the temperature and the pressure. They also quoted that the crystallization was inhibited for temperatures lower than T=573 K. Grzechnik, Grande and Slolen [35] studied the inverse process, the pressure induced amorphization of the GeSe 2 with Raman scattering and energy-dispersive X-ray diffraction. Their experimental work showed that no phase transition to known polymorphs of GeSe 2 was observed until 11GPa. However, between 11 and 14 GPa, the sample showed a darkening related to a possible transition from 2:4 to 3:6 coordination during compression. At room temperature, the transition from four fold to six fold coordination is kinetically inhibited in crystalline materials, and in the range of 7-14 GPa a large disordered phase forms, the β -GeSe 2 is no longer stable. It was concluded that the pressure induced amorphization of AB 2 compounds results from a tendency to reach a higher coordinated phases at different pressures depending on the parent structure as well as the ionicity of the system. More light was brought to the GeSe 2 crystalline and glassy phases by Grande and Co. [36] using X-ray diffraction measurements under pressure up to 7.7 GPa and temperature up to T = 1100 • C. GeSe 2 was found to be in a 2D crystalline structure where the ambient pressure structure is made of 50% of edge and 50% of corner sharing Ge. In contrast to Shimada work, it was found that amorphous GeSe 2 compressed at T = 250 • C and T = 300 • C was observed to crystallize at 3.0 and 7.7 GPa. Also, during quenching from 7.7 GPa and 500 • C an orthorhombic polymorph of GeSe 2 was found.

Further study by the Grande team [37] focused on the crystalline GeSe 2 by using angledispersive X-ray diffraction at the pressure range 2-6 GPa. This study confirmed the transition to a three-dimensional crystalline structure above 2GPa and at 698K. This transition is ex-plained by anisotropic lattice distortion due to the cooperative tilting of rigid GeSe 4 tetrahedra. The authors claimed that a possible transition in the glassy framework of the germanium dichalcogenides at high pressure and high temperature (including an anomalous compressibility and thermal expansion) could be similar to the one in the crystalline GeSe 2 . In analogy to crystalline polyamorphism, this behaviour could be called "amorphous polymorphism".

In addition to the crystalline germanium diselenide, the liquid state had been also subject to X-ray diffraction study under pressure [38]. It was found that under compression the connectivity of the liquid changes from a 2D to 3D. The 2D structure is characterized by a mixture of edge and corner sharing tetrahedra, while the 3D one is only made by corner sharing tetrahedra. This conversion from a mixture ES/CS to a full CS network is accompanied by breakdown of the IRO, completely lost above 2.5 GPa.

The occurrence of a structural transition on the liquid and the crystalline states, in addition to the two amorphous structures which mimic the 2D to 3D transition observed by Shimada, suggest that a first order transition between two liquid GeSe 2 with different densities could be possible. Furthermore, measurements based on high pressure Raman spectroscopy on the glassy GeSe 2 [39] showed that the Raman spectra are reversible. In addition and up to 9 GPa, no structural transition was reported and only a change in the sample color was noticed. The topological changes in glassy GeSe 2 has been investigated by high-energy x-ray and neutron diffraction at pressures up to 9.3 GPa [40] (e.g see Fig. [2.12]).

These measurements interpreted the densification mechanism as a continuous increase of the mean coordination number with a conversion from edge to corner sharing tetrahedra. The diffraction data showed that the conversion is accompanied by an elongation of the mean Ge-Se bond length and a substantial distortion of the GeSe 4 tetrahedra. Moreover, by following the changes on the pair correlation functions as a function of pressure, the IRO showed a decrease with increasing the pressure to 9.3 GPa accompanied by a decrease in the intensity of the FSDP. Recently, glassy GeSe 2 has been subject of acoustic measurements using synchrotron radiation under pressure up to 9.6 GPa [42]. The outcome of these measurements showed a minimum in the network rigidity around 4GPa. The system responds differently depending on the applied pressure, for pressures until 3GPa the ES convert to CS maintaining the flexible character of the network. However, for higher pressure, a continuous increase on the density and thus the coordination persists, while the ratio ES to CS remain constant Fig. [2.13].

As summarized above, many experimental efforts were devoted to study and understand the behaviour of disordered phases of GeSe 2 under pressure. From the theoretical point of view the efforts deployed so far to understand the behavior of these glasses under pressure are by far not comparable to the availability of experimental data. To date, the only work on this area is the one by the team of Drabold [41]. In that work, the structural changes of the amorphous GeSe 2 under pressure up to 75GPa were investigated. The simulation were performed by using the Parinello-Rahman method and within the framework of the non self-consistent first-principles molecular dynamics. This method employs DFT within the local density approximation and the Harris functional. In agreement with the experimental work, this study revealed a continuous and reversible structural transition upon compression into a metallic amorphous phase where the Ge and Se coordination's change gradually from 2 and 4 to 4 and 6, respectively. For high pressures, Se clusters appear in a chain like fashion. The FSDP was found to shifts to large k values with a decreasing intensity. Regarding the structure, the mean bond length does not change considerably until 12GPa while in the range of 12-16 GPa it increases substantially.

Similarly, the mean coordination number increase until 12GPa. The pressure was found to suppress the homopolar Ge-Ge bonds.

The nature of the density-driven network collapse in GeSe 2 glass is elusive. The different set of results at pressures up to 9 GPa could be interpreted in terms of (i) a continuous structural change in which the mean coordination number n steadily increases [40,41] as CS motifs are initially replaced by ES motifs [1], or (ii) a two steps mechanism [42], before 3 GPa the network is found to be flexible as the ES are replaced by CS motifs, beyond the threshold of 3GPa the mean coordination number n increases leading to a rigid network, or more simple (iii) the network topology remains unchanged [39].

As in the case of glassy GeSe 2 , the structural changes under pressure of amorphous GeSe 4 had been tracked by high-energy x-ray diffraction experiments up 8.6 GPa [43]. The main outcome of this work rely on finding an invariant mean coordination number with increasing the pressure. This result is in striking contrast with the case of GeSe 2 where a gradual increase of the coordination was found (e.g see Fig. [40]. The blue triangles are ambient data from Ref. [44]. Figure from Ref. [43] The densification mechanism in amorphous GeSe 4 was interpreted in terms of the shifts of the second neighbor and the higher coordination shells. The peaks related to these shells in the pair correlation functions move to shorter distances and increase in intensity indicating a tendency of the network to fill the voids present at ambient pressure. Within the Ge x Se 1-x family, the origin of the densification mechanism depends on the bonding and the chemistry of the network. Many possible interpretations were proposed to explain and describe the transition mechanism. In order to discriminate between all the proposed mechanisms and provide a clear atomic scale description of the related structural and electronic changes under the pressure effect, it highly desirable to perform first-principles self-consistent simulations. Given the above context, in the present thesis work, two prototypical binary glasses within the family Ge x Se 1-x have been selected as target systems to perform intensive first-principles molecular dynamics simulations under pressure: glassy GeSe 2 and glassy GeSe 4 . amorphous systems as well as all the related details are detailed in appendix A.

Enhancing ambient pressure model

Before looking into the changes in the GeSe 2 network topology under pressure, it is of great interest to make sure that available models (studied at fixed volume, corresponding to the experimental density) are indeed available at a vanishing (ambient) pressure. Previous theoretical studies carried out at the experimental density, under the assumption of ambient pressure, were instrumental in providing a clear picture of the network structure. Briefly, the structure consists of a variety of structural units, encompassing GeSe 4 tetrahedra, homopolar bonds, and defective Ge and Se coordinations, accounting for a moderate departure from chemical order. Analysis of the methodology employed within FPMD reveals that most simulations have been performed at a fixed density. Even though these simulations are intended to be stress-free (zero pressure) no information on the actual values of the pressure were available. Account of the pressure is worthwhile since the density vs pressure relationship holding experimentally might not be exactly reproduced by the DFT model [2]. As a consequence, undesirable pressure effects altering the room temperature equilibrium properties can be observed. The starting point of this chapter corresponds to the observation that a set of FPMD simulations were carried out at 300K for a number density ρ a = 0.034 Å -3 , slightly larger than the value quoted in the experimental work ρ exp = 0.0334 Å -3 [3, 4]. In order to evaluate the pressure on the cell we have simulated a new system at the experimental density ρ exp , in addition we have also computed the pressure on the previously simulated system at ρ a . A residual pressure of ∼ 1GPa has been found in both model, showing that the model was not entirely stress-free for a volume (density) equal to the experimental value, this latter corresponding to ambient pressure. By expanding the simulation box and further relaxing the whole structure, we were able to lower the pressure to less than ∼ 0.1 GPa. The new corresponding density is ρ b = 0.0326 Å -3 . The final structure at the relaxed density was selected to start a thermal annealing cycle as described in the Annex A. Releasing this residual pressure has an important effect on the structural rearrangement. Figure [3.1] shows the new set of partial pair correlation functions compared to previous simulation results, simulation at the experimental density and the experiment results. The experimental results of Ref. [5] (solid curve) are compared with the results of Ref. [1] with ρ a (dashed red curve) and with the new simulation at the experimental density ρ exp (dotted blue curve with stars) and the present relaxed structure results at ρ b (NVT ensemble, dashed-dotted green curve).

Regarding to the previous simulations and experiment, the present g αβ for both ρ b and ρ exp present a small shifts in the peaks positions within the typical statistical fluctuations of the order 0.03-0.05 Å (see table [3.1]). When looking at the general shape, g GeSe and g SeSe are less affected by the change of the density, while g GeGe exhibits large changes. Adjusting the density to the experimental one enhances the shape of the third peak, and a clear minimum appears. The intensities of the first and the second peaks remain always largely over-estimated. The relaxed density model (ρ b ) provides a better agreement with the experimental g GeGe . The discrepancies between the second and the third peak intensities and the measured are considerably reduced, since the second peak intensity decreases and the third intensity peak increases. A quantitative comparison between these peaks could be provided by the partial coordination numbers n αβ , as obtained by integrating the corresponding g αβ peaks (see table [3.1]), also a striking improvement of these numbers is found for ρ b . Results for previous simulations at ρ a from Ref. [1] and experiment from Ref. [5] compared to new simulation at ρ exp and at the relaxed structure ρ b .

g αβ FPP (Å) n αβ IR (Å) SPP (Å) n ′ αβ IR (Å) TPP (Å) n " αβ IR (Å) g exp GeGe (

Amorphous germanium diselenide under pressure

Not only the simulation at the experimental density provide a better agreement with the measurements than the old ρ a one, but also the different coordination numbers as obtained from ρ b simulation feature provides an unprecedented improvement. The three main peaks in the g GeGe refers as r increases to Ge atoms involved in homopolar bonds, edge-sharing and corner-sharing connection respectively. The following table present the different Ge environment fractions as obtained from direct analysis:

N Ge-Ge (%) N Se-Se (%) N Ge (ES) (%) N Ge (CS) (%)
Simulation at ρ a Ref. [ Note that in Ref. [6], a molecular dynamics approach was used in conjunction with a reverse Monte Carlo method.

The current results within the relaxed model (ρ b ) shows a very good reproduction of the experimental fraction of Ge atoms involved in homopolar bonds. The fraction of corner sharing Ge atoms is improved upon previous simulation, this stems from the substantial change of the shape of the third g GeGe peak.

The reciprocal space properties also are improved by releasing the residual pressure. Fig. [3.2] shows the Bhatia-Thornton concentration-concentration S cc (k) partial structure factors. The correct reproduction of S cc (k) proved to be a stringent test of the FPMD approaches. Comparison between all the simulated models and the measurement shows that the calibrated system gives a better account of the S cc (k). The intensity of the main peak located at k ∼2 Å -1 is very close to the experimental one. This same height was underestimated when working at the ρ a . In addition a minimum appears in the range 4 Å -1 < k < 5.5 Å -1 , when simulated with ρ b . The new set of data improves the position and the shape of the first sharp diffraction peak (FSDP). for glassy GeSe 2 . Experimental results of Ref. [5] (black solid line), simulation results with ρ a Ref. [1](red dashed curve), simulation results with ρ exp (blue dotted curve and simulation results of the calibrated system at ρ b (green dashed curve). 91 3.1. Amorphous germanium diselenide under pressure Overall, the new calibrated system shows an enhanced chemical order with as many as 92 % of Ge atoms being fourfold coordinated, 72 % of them within a GeSe 4 tetrahedron. In the previous simulation, these values were equal to 78.1% and 62.5% respectively. In conclusion, account of residual pressure effects at room temperature leads to a better agreement between atomic-scale models and experiments. Therefore, the relaxed configuration at ρ b has been used as initial model to generate, study and predict the structural changes under pressure.

Topological changes under pressure

Equation of state

In order to cover the range of pressures reached in previous experiments and to go beyond, thirteen independent simulations were performed in the NVT ensemble and five independent simulations were performed in the NPT ensemble for pressures going up to ∼ 16 GPa. Our NVT data were compared to neutron diffraction results from Prof. Salmon's team. Experimentally, neutron diffraction with isotope substitution method (NDIS) was employed to investigate the structural organisation via the study of the partial structure factor. This latter technique has been applied using a D4C [7] diffractometer on 70 Ge N Se 2 and 73 Ge 76 Se 2 up to a 8.2GPa pressure. For higher pressures the structure of N Ge N Se 2 has been investigated using the diffractometer PEARL [8], (see appendix B for hints on the diffraction theory). In what follows and for the sake of consistency, simulation and diffraction data will be compared at equal reduced number density. Table [3.3] illustrates the different set of comparisons.

In figure [3.3], we plot the pressure dependence of the reduced volume. The computed equation of state from the NVT simulations is in very good agreement with X-ray diffraction results from Mei and coworkers [9] and MD simulations by Durandurdu [10], account being made for the error bars. The equation of state shows some discrepancy with the experimental one, since for the same pressure values the density found from the NPT results are underestimated. For this reason, caution must be exercised when comparing NPT results with experiments. Nevertheless, NPT results will be shown and discussed for the sake of completeness when needed.

NVT Simulation

D4C Experiment PEARL Experiment P(GPa) ρ sim (Å -3 ) ρ sim /ρ 0 P(GPa) ρ exp (Å -3 ) ρ exp /ρ 0 P(GPa) ρ exp (Å -3 ) ρ exp /ρ 0 0 0. The NVT simulation and experimental data were fitted by a second order Birch-Murnaghan equation of state of type:

P = 3B 0 2 V V 0 -7 3 - V V 0 - 5 3 (3.1)
B 0 gives the bulk modulus value that expresses the resistance of the compressed material to the applied pressure. The computed values from simulation and experiment are respectively; 11.93 GPa and 10.65 GPa.

Reciprocal space properties

Recently, progress on instrumentation and methodology of neutron diffraction allowed experimentalist to perform accurate measurements of the neutron diffraction patterns with isotopic substitution (NDIS) for multicomponent glasses and liquids at high pressures. This technique provides the different partial structure factors and the different difference functions obtained by computing the difference between the structure factors obtained from different isotopes (see appendix B for details). The partial structure factors for the different isotopes can be compared to those computed from the simulation. Hints on the diffraction theory and the computational details used to get the isotopic structure factors from theory and experiments are provided in appendix B. The agreement between the two sets of data is very good for all pressures. The main peaks are reproduced (intensities and positions) within the statistical errors. The statistical noise in the computed structure factors is due to the finite length of the analysed trajectory, a possible way to reduce this noise is to perform longer simulations or several, distinct, independent runs. One should note that these data are the "direct outcome" data of the reciprocal space calculation with no application of any smoothing or window function.

The computed structure factors show systematic changes with increasing the pressure. The FSDP is associated with oscillations of periodicity 2π/k FSDP in real space whose extent is controlled by the correlation length 2π/∆k FSDP , where ∆k FSDP is the full width at half maximum of the FSDP [11]. In our case, k 0GPa FSDP = 1.15 Å -1 and ∆k 0GPa FSDP = 0.55 Å -1 which gives a mean periodicity of ∼ 5.5 Å and a correlation length of ∼ 11.5 Å. The FSDP moves rightward and decreases in intensity when increasing the pressure, reflecting a significant change due to the onset of the intermediate range order. The main peak at ∼ 2 Å -1 at 0 GPa is associated with the extended range order [12]. This peak shifts to higher wave vector values and sharpens, indicating also a changes in the extended range order. For pressures until 8.2 GPa, the FSDP decreases in magnitude and disappears completely for the highest pressure. PEARL experiment at pressures ranging from 8.7 GPa to 16.1 GPa were not able to access the region of k < 1.55 Å -1 , the data plotted representing the fitted Lorentzian functions. Looking at the FPMD results, a vanishing bump persists and indicates that the IRO does not exhibit important changes for values higher than 8.5 GPa.

A further description is made available by the difference functions ∆F Ge (k) and ∆F Se (k) shown in Fig. [3.5]. Experimental results are only provided for pressures up to 8.2 GPa and compared to the ones constructed from our calculations in the same experimental fashion. The pronounced FSDP in the ∆F Ge (k) and its absence in the case of ∆F Se (k) confirms that this peak is mainly due to Ge-Ge correlations as already suggested in the literature. The FSDP magnitude follows the same trend seen in the total structure factor. However, unlike in the ∆F Se (k) case, the ∆F Ge (k) shows oscillations for wave vectors as high as ∼ 12 Å -1 for all pressures, reflecting the high structural order of the Ge sub-network. A steady decrease of the depth of the wall around 2 Å -1 is associated to possible changes altering the intermediate and short range order.

-0,5 0 0,5

Real space properties

The real space properties could be assessed through the study of the pair distribution functions. Two methods can be used. The first amounts to Fourier transforming the obtained structure factors from reciprocal to real space. In the second method, one calculates the partial or the total PCF's directly in the real space. By following the first method, the corresponding Fourier transform real space difference functions ∆G Ge (Q) and ∆G Se (Q) are shown in fig. [3.6]. A maximum cutoff k max = 19.55 Å -1 was applied as employed when obtaining the high-pressure neutron diffraction data.

The G(r) functions from experiment (solid dark (black) curves) were obtained by Fourier transforming the measured F (k) functions shown in Fig.[3.4] via the application of a Lorch modification function. The broken dark (black) curves show the extent of the unphysical oscillations at r-values smaller than the distance of closest approach between the centers of two atoms. The solid light (green) curves

show the Fourier transforms of the simulated functions shown in Fig. [3.4] after applying the same maximum cutoff value k max = 19.55 Å -1 as used for the neutron diffraction data.

A very good reproduction of the intensities and the positions of the different peaks over all the real space interval is found for pressures up to 8.2 GPa. Unphysical oscillations on the short range scale are due to the Fourier transform finite integration window. Surprisingly, for the higher pressures, the real space difference function has some discrepancies with respect to the measured one. The intensity of the first peak is much higher than the measured one and, for pressure beyond 10.9 GPa, the position of the second peak is shifted to higher r values. These discrepancies are likely to be due to the different thermal histories of the experimental and simulation samples. Indeed, in the diffraction experiments, glassy GeSe 2 was cold-compressed. Experiments were carried out through a progressive increase of the load. At each pressure, diffraction patterns were collected for some hours at fixed temperatures T=300 K before raising up the load. As any other glassy structure, amorphous GeSe 2 is characterised by vanishing diffusion at T=300 K, leading to a strong memory effect on the experimental samples when the pressure is increased. This means that, for each increasing value of the pressure, the structure of the system keeps some memory of the previous configurations at a lower pressure. On the other hand, in the modelling case (our work), we adopted a different strategy. We generated independent systems under pressure starting from the same initial point; each system undertook a thermal annealing cycle under pressure (see appendix A). Unlike the experimental situation, much higher temperatures are attained, characterized by substantial diffusion. In this way, our systems can explore the configuration space at high temperature under pressure before relaxing at room temperature. For the time being, no experiments featuring a thermal annealing cycle under pressure (as it was done within our simulations) can be performed. As a working hypothesis, the experimental and the modelling procedures can be reconciled by selecting the first part of the annealing run at T=300 K as the part of the full simulation cycle closer to the experimental cold compression procedure. Fig. [3.7] shows the comparison of the total pair correlation function obtained from ND experiment, simulations with annealing and simulation of direct compressed systems. Overall, the cold compression simulations results improve the agreement with the experimental results. The first peak position is either preserved or shifted to the smaller values. In addition, the first minimum positions are shifted rightward and for r > 5 Å. As a consequence, the discrepancies observed above while comparing the thermal cycle simulation results and the experiments could be associated to the different application of the pressure (either by introducing a thermal cycle or by directly compressing the sample at room temperature).

Results from direct calculations of the partial pair correlation functions (PPCF) are illustrated on Fig. [3.8]. As the pressure is increased, dramatic changes alter the different partial pair correlation functions. Two regions with two different behaviours could be distinguished. If we focus on g GeSe (r), a main peak around 2.35 Å is found, representing the typical Ge-Se bond length at ambient pressure. This peak remains at the same position and maintains its width up to 8.7GPa, with a small decrease in its intensity. For higher pressures, this same peak shifts rightward and broadens by decreasing in magnitude. For g GeGe (r), the three main peaks are 3.1. Amorphous germanium diselenide under pressure linked to homopolar, corner-sharing (CS) and edge-sharing (ES) configurations. The second peak intensity decreases dramatically for P > 8.7 GPa. The third peak moves leftward and broaden dramatically to include the second peak. For values larger than 8.7 GPa, only one broad peak persists, containing contributions from Ge involved in ES and CS conformations, now found at close Ge-Ge distances. Finally, g SeSe (r) exhibits also remarkable changes when the density increases. The first peak associated with the homopolar Se atoms broadens and shifts to higher distances, while at the same time the main peak representative of the second coordination shell shifts to a lower distance. This reflects that Se atoms in the GeSe 4 tetrahedra are closer to each other. The existence of two regimes in the structural behavior as a function of the applied pressure is related to a structural transition with increasing the density. In order to characterize the changes on the structural properties as the pressure is increased, local structural analyses are provided in the following paragraph.

Structural properties

Further insight into the network sensitivity to the pressure effects can be extracted from the local structure analysis. One could compute the different atomic bond lengths r αβ as given by the position of the first peak for each partials g αβ (r) where α and β are the chemical species. The partial coordination numbers nβ α of chemical species β around an atom of chemical species α is given by the integral of the corresponding g αβ (r) peak including distances up to the first minimum positions. The mean coordination number is given by n ≡ c Ge [ nGe Ge + nSe Ge ] + c Se [ nSe Se + nGe Se ] where c Ge and c Se are the concentration of Germanium and Selenium. From simulations, the n αβ can be derived by integrating the first peak of the corresponding g αβ (r) . Experimentally, n can be extracted from the first peak in the Bhatia-Thornton number-number partial pair-distribution function g NN (r), which can be accurately measured by making a neutron diffraction experiment on a sample of N Ge N Se 2 . In figure [3.9] we display a comparison between our results on the NVT and NPT ensembles and the experimental mean coordination number and the mean first peak position. The results of previous X-ray diffraction measurement [9] are also provided after removing an unphysical slope on the measured k-space functions and renormalizing.

samples on D4 and averaging the results [(green)]; (ii) x-ray diffraction [9] [(magenta)]; or (iii) FPMD NVT simulations [broken (black) curve], or cold compression NVT simulations at T= 300 K [dotted broken (maroon) curve with cross marks], or NPT simulations [dotted (dark green) curve with stars] or HSE06 simulation [orange cross mark]. In the bottom figure, the horizontal chained line gives the "8-N" rule expectation of n = 2.67.

For values smaller than P ∼ 8.5 GPa (ρ/ρ 0 ∼ 1.42), the mean nearest-neighbor peak position and the mean coordination number are in excellent agreement with the measurements within the experimental error bar. The mean coordination number is consistent with the 8-N rule and the network retains its ambient pressure topology. The fractions of the CS and ES Ge as shown in Fig. [3.10] exhibit small fluctuations as the pressure increase. At higher pressures, the measured and calculated results show an increase in n as the first coordination shell expands to incorporate a larger number of neighbors, although this process starts at 8.5 GPa from NVT results (ρ/ρ 0 ∼ 1.42) as compared to 12GPa (ρ/ρ 0 ∼ 1.55) from experiment. The simulated direct compression results at ambient temperature lead to a better agreement with the experiments. The mean coordination number is smaller than the one found from the simulation with the thermal cycle. This result is consistent with the proposed explanation of the discrepancies experiment/theory at high pressures. A computational approach using the HSE06 hybrid density functional [13] was employed with a single pressure point of ∼9 GPa in the aim of enhancing the XC description by retaining a more ionic and less metallic character for the chemical bonding [14]. Calculation details and some further results are provided in appendix E. A good agreement is found for the mean co- ordination number with the BLYP simulations. In contrast, the mean bond length is by far over estimated. The hybrid scheme fails to reduce the gap experiment/theory at high pressures. By looking at the NPT simulation results, the structural changes follow the same trends. However, the pressure threshold for the occurrence of a structural transformation cannot be determined because of the limited number of different applied pressures. Despite this limitation, the simulation r values are in a good agreement with the X-ray results, while they do not agree with the neutron diffraction results, exhibiting only a slight increase of the mean coordination number. This amounts to concluding that, in spite of any consideration on the intrinsic quality of our theoretical data, we are in the presence of two sets of experimental results essentially not in agreement with each other. As an additional comment, we note that the simulated direct compression r are in a good agreement with the NDIS data. The fraction of Ge CS reduces and give place to Ge ES as the pressure increase. The breakdown of the Ge ES into Ge 1 , Ge 2 and Ge 3&4 (see appendix C for definitions) highlights the progressive appearance of atoms of type Ge 2 followed by atoms of type Ge 3&4 . For the highest pressures, the network become fully connected in a edge-sharing fashion, where the Ge atoms reach a maximum coordination of six and the Se a maximum coordination of four. This supports the occurrence of a phase transition around 8.5 GPa from tetrahedral to octahedral network.

Local environment analysis

The pressure dependence of the intra-polyhedral Se-Ge-Se bond angle (θ SeGeSe ) and the interpolyhedral Ge -Se -Ge bond angle (θ GeSeGe )) distributions as computed from our simulation are shown in Fig. [3.11]. The bond angle distributions were computed as follows: The pressure dependence of the intra-polyhedral Se-Ge-Se and the interpolyhedral Ge-Se-Ge bond angle distributions. NVT simulations (left) and NPT simulations (right).

P(θ ) = B(θ )/sin(θ ) (3.2)
Two main peaks appear in the Ge -Se -Ge bond-angle distribution (BAD), the first one is ascribed to edge-sharing Ge centred motifs sharing two Se atoms, while the second one is attributed to corner-sharing Ge centred motifs which share only one Se atom. At NVT ambient pressure this two peaks occur respectively at 79˚and 103˚, in agreement with previous experiments by Salmon [12] where the peaks were found to take place around 80˚and 98˚. The Se -Ge -Se bond angle distribution features one main peak around 111˚, due to the predominant tetrahedral nature of the network. By increasing the pressure up to ∼ 8.5GPa, the BAD's are less affected and the peak position do not change. Beyond this threshold pressure, all the peaks position shifts leftward to reach 90˚and in the Se -Ge -Se BAD a new peak appear around 180˚. The picture given by NVT simulations has to be complemented by the one due to NPT simulation. The Ge -Se -Ge distributions conserve the peak positions and exhibit a moderate broadening with increasing pressure. Only for ∼ 16GPa the peak in the Ge -Se -Ge distribution merge and get centered around 90˚. The Se -Ge -Se peak shifts to 90˚simultaneously with the appearance of a small peak around 180˚. Fig. [3.12] shows the comparison between the NVT mean bond angles θ GeSeGe and θ SeGeSe 104 105 3.1. Amorphous germanium diselenide under pressure computed by the following expression with an upper integration limit of 140˚:

θ = dθ θ P(θ )/ dθ P(θ ) (3.3)
One could estimate the Se -Ge -Se and Ge -Se -Ge bond angles from the experiment by using the cosine rule:

cos(θ SeGeSe ) = 1 -r 2 SeSe /2r 2 GeSe (3.4) cos(θ GeSeGe ) = 1 -r 2 GeGe /2r 2 GeSe (3.5)
An accurate estimate of the nearest-neighbor distance is required to get a good estimate of the bond angles. The simulation tracks successfully the shift of the peak positions with increasing density (pressure) found experimentally. The offset between the simulated and the measured Ge-Se-Ge bond angles is due to the poor accuracy in measuring the Ge-Ge distances. The network is exhibiting a phase transition from perfectly tetrahedral arrangement to a new pseudo cubic one where the different αβα angles are centered around near 90˚. In order to elucidate the new atomistic This analysis is supported by the rings analysis (see Fig. [6.19]). Up to 9 GPa the network is dominated by 4-, 5-and 6-membered rings with a very small fraction of large rings. Beyond 9 GPa, the structure exhibits an increase of the 4 to 6 membered rings, consistent with the phase transition from tetrahedral to defective octahedral configurations. The counts of large rings increases also, but it remains very small compared to the amount of 4-fold rings. The density dependence of the fractions of n-fold coordinated Ge and Se atoms (n = 2, 3, 4, 5 or 6) is provided in Fig. [3.15], the breakdown of these fractions into the proportions of these n-fold species that contain homopolar bonds is also provided in the same figure (see appendix C for definitions). The evolution of these fractions show again the two phases behaviour, up to ∼ 8.5 GPa the 4-fold Ge and 2-fold Se dominate in agreement with the tetrahedral arrangement. Less than 10% of 3-fold Ge and 3-fold Se exist. In the second stage, the fraction of 4-fold Ge (2-fold Se) decreases progressively and is replaced by 5-and 6-fold Ge (3-, 4-and 5-fold Se). The occurrence of the high coordinated regime is in agreement with the increase of the mean The pressure dependence is given for the fractions of n-fold coordinated Ge and Se atoms, along with the fractions of these n-fold coordinated Ge and Se atoms that contain homopolar bonding defects. The error bars (usually smaller than the symbol size) were calculated according to Ref. [16]. Left plot: NVT simulation results, right plot: NPT simulation results. The NPT results show again the same trends found by the NVT simulations, however, for the highest reached pressure the fractions of the overcoordinated Ge and Se atoms are very small and the transition occurs for higher pressure.

The breakdown of these n-fold atoms into defected n-fold homopolar atom shows that at the first appearance of the overcoordinated atoms they contain many homopolar connections. The homopolar atoms, thus, play a crucial role when the transition is about to happen.

Structural transition mechanism

The density-driven network collapse of the structure could be described as a two stage transition. In the first stage, as the density is increased to ρ/ρ 0 ∼ 1.42 (P ∼ 8.5 GPa), the mean coordination number remain constant, the chemical order is maintained and the network topology is marginally affected. Nevertheless, small fluctuations of the number of CS and ES of tetrahedral units occur. The mean inter-tetrahedral Ge-Se-Ge bond angle decreases from ∼103t o ∼97˚for CS units but remains invariant for ES units at ∼79˚. These results are in contrast with the case of the crystalline phase, where transitions are observed from an ambient pressure 2-D structure (in which layers are formed from equal numbers of ES and CS tetrahedra), to 3-D structures of densely-packed CS tetrahedra. Also, this description is in contrast with the one proposed by Drabold [10] and Mei [9] where the densifcation mechanism was interpreted in terms of continuous increase of the mean coordination number. In a second stage, as the density is increased beyond ρ/ρ 0 ∼1.42, n and r increases steadily. In the latter, four-fold coordinated Ge atoms become 5-and 6-fold coordinated Ge atoms and 2fold coordinated Se atoms transform to higher coordinated Se atoms where homopolar bonds are common (see Fig. [3.15]), and there is a monotonic reduction in the CS to ES ratio (Fig. 3.10). The FPMD results indicate an increased metallic character with pressure, whereas the diffraction results for n are consistent with the retention of semi-conducting behavior. A computational approach using the HSE06 hybrid density functional [13] was therefore employed with the aim of retaining to higher pressures a more ionic and less metallic character for the chemical bonding [14]. The calculations for a single pressure point of ∼9 GPa did not, however, lead to improved agreement with experiment. The origin of the discrepancy may therefore lie in the exploration of different pathways of meta-stability in a pressure regime where the structure of GeSe 2 glass can rapidly change: the glass was cold-compressed in the diffraction experiments but was annealed at pressure in the FPMD simulations. Indeed, a discontinuous semiconductorglass to metal-crystalline transition has been reported for GeSe 2 at P ≃7 GPa [17,18]. (P ∼ 4 GPa) [19]. A minimum was inferred for the network rigidity at P ∼ 4 GPa, originating from a competition between two densification mechanisms: a conversion from ES to CS tetrahedra which initially enhances the network flexibility versus an increase in n which stiffens the network. An increase in the CS to ES ratio from 1.3 to 1.7 in the range ρ/ρ 0 ∼1.07-1.20 (P ∼1.29-3.55 GPa), consistent with the trend observed by Raman [20] spectroscopy is found from the FPMD results. However, since the present work shows that n does not increase until ρ/ρ 0 1.42 (P ∼ 8.5 GPa), the stiffening mechanism is most likely related to more subtle changes in the network connectivity e.g. to an adjustment of the CS to ES ratio at ρ/ρ 0 ∼1.24 (P ∼4 GPa). Fig. [3.16] shows typical atomistic configurations for different points in the densification process. When ρ/ρ 0 increases beyond 1.42 (∼ 8.5 GPa), the Se atoms in Ge 0 motifs (CS mo- tifs) make additional bonds with their neighboring Ge 1 motifs, transforming the latter to 5-fold coordinated Ge 2 units having a distorted square pyramidal geometry. As Ge 0 tetrahedra are removed, the Se-Ge-Se bond angle distribution broadens, its peak shifts from ∼111˚towards 110 111 3.1. Amorphous germanium diselenide under pressure smaller angles, and there is a merging of the twin peaks in the Ge-Se-Ge bond angle distribution (see figure [3.11]). Around ρ/ρ 0 ∼ 1.53 (15.2 GPa), the fraction of Ge 0 units continues to de- crease, leaving their place to more Ge 2 and a few Ge 3&4 type configurations, where many of the latter correspond to 6-fold coordinated Ge atoms. Accordingly, the main peak in the Se-Ge-Se bond angle distribution shifts to around 90˚and a second peak appears near 180˚. Eventually, the network develops a pseudo-cubic arrangement of Ge-centered units for which the main peak in the Ge-Se-Ge bond angle distribution is found around ∼90˚. In all of this, the fractions of defected 5-and 6-fold coordinated Ge atoms shows a maximum around the pressure for which they start forming. In addition the fractions of defected 4-and 5-fold Se coordinated atoms increase rapidly beyond the structural threshold pressure. These findings allow to conclude that the homopolar bonds mediate in the initial development of the higher-coordinated Ge-centered polyhedra.

Amorphous GeSe 4 under pressure

The initial configuration of GeSe 4 was taken from ref [21] and the simulation details are exactly the same as those used to simulate the GeSe 2 .

The GeSe 4 glass belongs to the Se rich side of Ge x Se 1-x family, since it corresponds to x = 0.2. For this composition, several concepts related to the phenomenological model of constraints can be employed. Accordingly, at this composition, the network exhibits a change from an elastically floppy or under-constrained network to a stressed-rigid or over-constrained network [20,22].

Equation of state

The simulations were performed both in the NPT and the NVT ensembles. Five simulation points covering pressures ranging from ambient to 13.2 GPa were generated. Our results will be compared to X-ray experiment recently published by Skinner [23] and, when needed, to Kalkan's results [24]. The values for which experiments and simulation will be compared are presented in table [3.4]

Reciprocal space properties

The dependence on pressure of the measured and simulated total X-ray structure factor is presented in Fig. [3.18-b]. A comparison of the first sharp diffraction peak position (FSDP) and the principle peak position (PPP) is also given in Fig. [3.18-a]. Our simulations reproduce the FSDP pressure dependence, with a non-negligible statistical noise found in the region 1.5 Å -1 < k < 4 Å -1 . For this reason, the Fourier transform of the real space pair correlation function is also provided in Fig. [3.18-b]. The Fourier transformed data show a non physical contributions in the low k region. At ambient pressure, the FSDP occurs at k = 1.13 Å -1 in agreement with X-ray measurements. Compared to the GeSe 2 case, the FSDP occurs at the same position (k = 1.15 Å -1 ) with a lower intensity and a smaller width. These differences could be due to the smaller number of Ge atoms in GeSe 4 , since the FSDP in the case of GeSe 2 was attributed to contributions of real-space Ge-Ge correlations over the intermediate length scale. For k ≥ 2Å -1 the peak positions and intensities are correctly reproduced from the Fourier transformed S X-ray (k). Figure [3.18-a] shows the pressure evolution of the FSDP and the principle peaks positions. As the pressure is increased, both peaks shift to higher wave vectors values and slightly sharpen. The principle peak intensity increases steadily and only for the highest simulated pressure it decreases. The changes in the FSDP intensity and position indicate a significant changes in the IRO with increasing the pressure. More insights into the reciprocal space properties are obtained through the Faber-Ziman partial structure factors. The data presented on figure [3.19] The pressure evolution of the partial Faber-Ziman structure factor from from Fourier transform of the partial real space counterpart presented on figure. [3.20]. From the top to the bottom S FZ SeSe (k), S FZ GeSe (k) and S FZ SeSe (k).

from Ge-Ge and Ge-Se correlations as proposed in th past [5]. As the pressure is increased, the FSDP in F FZ GeGe (k) decreases and shifts rapidly to higher k values, while the changes are less abrupt in F FZ GeSe (k). The principle peak located around k ∼ 2 Å -1 in the total structure factor results from Ge-Ge and Se-Se correlations, somewhat diminished by a negative Ge-Se contribution. When moving to the high pressure regime, the principle peak intensity as found in the total S X (k) increase, this increase being mainly due to the F FZ SeSe (k) principle peak intensity.

Real space properties

The pressure evolution of the total and the different partial pair distribution functions are shown in Fig. [3.20]. The main peak in the total g(r) occurs at 2.36 Å, in excellent agreement with experimental results [24]. This peak is due to Ge-Se and Se-Se correlations. With increasing the pressure, the first peak in g SeSe (r) and g GeSe (r) shift slightly to lower r values and causes the 2.9% shift of the PPP in g tot (r) at 13.2 GPa. This shift is comparable to the shift found in experiments (1.2% for 15.3 GPa). It was noticed from experiments that there is a small elongation of the Se-Se bonds. In addition, the measured peak intensity decreases, leading to a reduction in the intensity of the main peak in g tot (r) by 55%. The second peak in the total g(r) is located at ambient pressure at 3.82 Å (3.73 Å from experiment at 1.3GPa). This peak is a combined contribution of Ge-Ge and Se-Se correlations. As the pressure is increased, this peak shifts to 3.52 Å and a small shoulder appears near 3.2 Å (3.23 Å from experiment). This shoulder is due to the appearance of a peak on the g SeSe (r) at the same position. The appearance of this peak in the high density regime has to be ascribed to the deformation of Se-Se-Se chains upon compression, leading to a rearrangement of the second shell of coordination around the Se atoms. In the case of g GeGe (r) at ambient pressure, two dominant peaks exist around 2.97 Å and 3.68 Å respectively. These peaks are due to Ge atoms belonging to edge or corner sharing configurations, respectively. In contrast to GeSe 2 , we notice the absence of the first peak associated to homopolar Ge bonds. As the pressure increases to ∼ 8GPa, the first peak shifts slightly to 3.02 Å, while the second one broadens and shifts to 3.44 Å. For higher pressure, these two peaks merge to give one broad peak centered around 3.42 Å. For 13.2 GPa some homopolar motifs start to appear. A double peak at 2.35 Å and 2.58 Å is a signature of their existence. Since the nearest neighbor distances exhibit small changes over the range of compressions, the densification mechanism is substantially different from the one holding for glassy GeSe 2 .

In the present case, the densification occurs through a reorganisation of the second and higher coordination shells around the Se atoms. This hypothesis is consistent with the experimental findings [23,24].

Structural and local environment properties

Figure 3.21 Top: Ge atoms bonding to homopolar (black circles), corner-sharing (left green triangle) and edge sharing (blue right triangle) configuration. Se atoms belonging to homopolar connections (red square). Bottom: coordination numbers of Ge atoms (red triangles) and Se atoms (green triangle). Mean coordination number from FPMD (blue diamond) and experiments (black stars).

Fig [3.21] shows the pressure evolution of the fractions of Ge ES, CS, homopolar and Se homopolar. The fraction of homopolar Se bonds reduces slightly upon compression, most likely because of the breakdown of the long Se chains when higher pressures are reached. The fractions of Ge ES and CS show a small reduction, these small fluctuations are consistent with the change in the first and second peaks in the g GeGe (r). The tetrahedral arrangement is thus preserved and only the inter-tetrahedral connections are affected by the compression.

The partial coordination numbers computed by integrating the first peak in the related g(r) remain constant for all pressures. The resulting mean coordination number remains constant, in excellent agreement with Skinner's results [23]. Our results show a very small increase of this value from 2.38 at ambient pressure to 2.48 at 13.2 GPa, to be compared to an increase from 2.35 at ambient to 2.45 at 8.6 GPa from X-ray experiments. Both results are consistent with an increase of less than 5%, in contrast to the Kalkan's measurements [24] pointing to an increase from 2.5 at 1.3GPa to 2.7 at 10.8 GPa (increases by ∼ 8% ). The difference between the later measurement and our results could be assigned to the low k max values used to integrate Kalkan's 118 3.2. Amorphous GeSe 4 under pressure structure factors. Further details on the connectivity profile are shown in Fig. [3.23]. The stability of the structure with increasing the pressure manifests itself through the conservation of tetrahedral arrangements, the fraction of 4-fold Ge and 2-fold Se remaining dominant. For the highest pressure, these fractions decrease slightly and a small fraction of 5-fold Ge and 3-fold Se atoms do appear. The fraction of 2-fold Se are broken down into contribution of Ge -Se -Se, Ge -Se -Ge and Se -Se -Se motifs, referred to as AB, AA and BB units. Figure [3.22] exemplifies each one of these configurations. At ambient pressure, the fraction of AB units is predominant (∼ 44 % ), while the proportions of AA and AB units are slightly smaller with ∼ 28 % and ∼ 27 % respectively. This results is in agreement with our previous FPMD calculations in the NVT ensemble [21] and reflect the network connectivity and the absence of a phase segregation. The high fraction of homopolar Se bonds and the existence of BB motifs are consistent with the description of GeSe 4 tetrahedra linked by Se n bridges. As the pressure increases, the fraction of AB units remains dominant but it decreases considerably (23% for 13.2 GPa). Instead, the fractions of the AA and BB motifs decrease slightly indicating a breakdown of long Se chains.

In order to complete the atomic scale description of the topological changes under pressure of the GeSe 4 glass we present in figure [6.7] the pressure dependence of the A-B-A bond angles distributions. The Ge -Se -Ge bond angle distribution has two peaks located at 78.5˚and 102.5˚associated to the Ge atoms belonging to edge and corner sharing conformations respectively. As the pressure is increased these two peaks merge. For 13.2 GPa, the distribution takes the form of a broad peak centered around 92˚. The Se -Ge -Se distribution is characterized by a main peak located around 111˚. This distribution broadens and shifts to low θ values (107˚) for highest pressure. The Ge -Se -Se BAD shows small changes when increasing the pressure. The main peak shifts by ∼4˚from ambient pressure to 13.2 GPa. However, this distribution broadens, reflecting the change in the AB fraction and local environments. Our Se -Ge -Se and Ge -Se -Se BADs are different than the one obtained from EPSR [24] where a features were found around 60 ˚. The crucial role played by the Se chains is reflected in the Se -Se -Se BAD were the distribution has a peak moving from 107˚at ambient pressure to 96˚at 13.2 GPa. This change indicates the reduction of the angles formed at ambient pressure, corresponding to a twist of the Se chains in the attempt to fill the intertetrahedral space. This shift in the Se -Se -Se distribution is consistent with the split of the second peak of the g SeSe (r) at high pressure. The bump of the BAD Se -Se -Se found around 180˚as the pressure increase is correlated to the increase of the fraction of 3-fold Se. As the pressure increases, the Se atoms make additional bonds with other Se atoms and the 3-fold centered Se atoms are found to be in a defective tetrahedral like environment with Se -Se -Se angles near 96˚and 180˚.

In summary, the behaviour of the GeSe 4 under pressure differs substantially from the one of small changes as increasing the pressure. The tetrahedral network is preserved. However, the densification occurs through changes of the second and higher shells of coordinations around Se atoms, in line with what was shown in previous experimental studies. Figure [3.25] illustrates a comparison between the ambient and the high pressure phases. Pressure induces the formation of Ge-Ge homopolar bonds (green arrows), Se 3 (blue arrow) and Ge 5 (red arrow) for 13.2 GPa.

αβα configurations tend to form right angles (90˚) and chain like motifs (180 ˚) reflecting the strong change from the perfect ambient pressure GeSe 4 tetrahedra. Se-Se-Se chains bends thereby provoking a head-tail contact effect between Se atoms, corresponding to the appearance of a new peak in the g SeSe (r) at 3.2 Å.

Chapter 4. Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4 126 presented density functional theory as a powerful tool to solve electronic structure problems in condensed mater. However, the structures of Te-based materials, as obtained by diffraction experiments, are not always accurately reproduced within DFT methods [1][2][3][4][5][6]. The major drawbacks of DFT regarding the Te-based alloys are the overestimation of the atomic bond lengths and the inaccuracies in the description of the atomic bonding. The account of the van der Waals (vdW) dispersion forces has been found promising for oxides [7], solids [8,9], water [10] and very recently in the case of the Ge 15 Te 85 liquid structure [11]. In this latter work, the account of the dispersion forces resulted in a substantial improvement of the atomic bond lengths and in an increased level of agreement between the computed and measured structure factors and radial correlation functions. Along the same lines, I have investigated the effect of the dispersion forces on the structural properties of amorphous GeTe 4 .

Simulation details

Amorphous GeTe 4 was generated (ρ=0.030186 Å -3 ) and studied by mean of first principles molecular dynamics as implemented in CPMD package. The computational scheme is the same one used for studying the Ge x Se 1-x family, with a larger plane wave expansion energy cutoff (here E c = 40Ry) and an integration time of ∆t =0.17 fs. The empirical vdW corrections was implemented as proposed by Grimme [12]:

E disp = -s 6 N at -1 ∑ i=1 N at ∑ j=i+1 C i j 6 R 6 i j f dmp (R i j ) (4.1)
Where N at is the number of atoms in the system, C i j 6 and R i j are the dispersion coefficient and the interatomic distance for an atomic pair i j. s 6 is a global scaling factor depending on the adopted density functional scheme. f dmp (R i j ) is given by:

f dmp (R i j ) = 1 1 + e -d(R i j /R r -1) (4.2)
This damping function is introduced to avoid singularities at shorter distances. In our case we used d = 20 and a cutoff distances R GeGe = 3.45 Å, R GeTe = 3.62 Å and R TeTe = 3.78 Å. The dispersion coefficients are given by C GeGe = 176.81 eV.Å 6 , C GeTe = 240.89 eV.Å 6 and C TeTe = 328.19 eV.Å 6 . These coefficients are calculated as C i j = C i C j where C i are the dispersion coefficients according to the London formula [12].

Impact of the dispersion forces

By taking as initial configuration the one obtained by changing the chemical identities of a GeSe 4 amorphous model, two different systems were produced. The first model consists of a 215 atoms simulated with the BLYP exchange and correlation functional by accounting (called hereafter BLYP+vdW) or not accounting (called hereafter BLYP) for the vdW corrections. Both systems (BLYP and BLYP+vdW) were subjected to a thermal annealing cycle (BLYP+vdW: 4ps @100K, 12ps @300K, 22ps @600K, 17ps @900K, 20ps @600K and 31ps @300K; BLYP: 2ps @300K, 12ps @300K, 19ps @600K, 20ps @900K, 20ps @600K and 32ps @300K). To consider also another exchange-correlation functional [13][14][15], we performed a second simulation on a 185 atoms model using the PBE XC functional with (called hereafter PBE+vdW) and without (called hereafter PBE) empirical corrections. These calculations were performed by a collaborator in Rennes (Dr. E. Furet) by following a procedure similar to the one described above.

Impact of the dispersion forces

The total structure factor and total pair correlation function are displayed in Fig. [4.1] and compared to the experimental results from ref. [5]. If one focuses, for the entire k range, on the computed S(k) without vdW corrections the agreement is fair in the PBE case and very good in the BLYP case. In particular, the first peak position is displaced rightward for PBE while, in the BLYP case, it is at the right position, with a good agreement for the intensity. After account of the vdW term, for both BLYP and PBE models, the first peak position get closer to the measured one and its intensity increases. For r > 4 Å -1 , the computed S(k) with and without vdW corrections exhibits a very similar patterns.

The impact of the vdW corrections on the real space properties could be ascertained by looking at the total pair correlation function. A priory, since the dispersion forces implemented by Grimme act basically on the range of r beyond R c , one would not expect a noticeable impact for short distances. Regardless of the choice made for the XC functional, for r > 3Å, a higher level of structural organisation manifests itself through the appearance of a deeper first minimum. In addition to a higher intensity, the second peak moves leftward.

It is worth to notice that only the model BLYP+vdW reproduces the experimental feature around 7.7 Å in the total experimental g(r) (knowing that most likely the patterns shown in the inset of figure [4.1] exhibit differences that might be of the order of the statistical errors). Overall, when looking at the two above properties in reciprocal and real space, the following conclusion can be drawn. In reciprocal space, vdW corrections enhance the intensities of the Chapter 4. Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4 128 peaks, with small changes for the peak positions. A moderate improvement can be found in the depth of the first minimum, even though the profiles features an additional bump absent when the vdW are not taken into account. In real space, the most striking change attributable to the vdW correction can be found in the shape of the total pair correlation function around the first minimum, definitely more pronounced. The peak positions and the partial coordination numbers, as computed from the partial pair correlation functions are displayed in the following Bond distances r i j (position of the first maximum of the pair distribution functions g i j ), partial coordination numbers n i j (integration of g i j up to the first minimum ), coordination numbers n i and the mean coordination number n.

For both the BLYP and PBE cases, g GeGe (r) exhibits appreciable changes with account of vdW corrections. In the BLYP+vdW case, the intensities of the first and the second peak increase. In the PBE+vdW case, the first peak is shifted toward the small r values and narrows. Also, a well defined minimum exists between the first and the second peak, corresponding to a well defined bond distance. These changes in the shape and position of the first peak induce changes on the n GeGe coordination number. As displayed in table [6.1], n GeGe for PBE+vdW and BLYP+vdW show opposite variation when turning on the vdW corrections.

The shape of g GeTe (r) does not change significantly after inclusion of the dispersion forces. Nevertheless, a small shift of the first minimum position leads to reduction of the n GeTe coordination number (see table [6.1]). The second peak intensities are slightly boosted for both models (BLYP+vdW and PBE+vdW). The small effects of vdW corrections on the Ge -Te bonds are ascribed to the strong iono-covalent nature of this specific bond.

On the contrary, among the three partial correlation functions, g TeTe (r) is the one undergoing the most noticeable changes when turning on the vdW forces. The Te -Te bond length reduces from 2.89 Å to 2.84 Å for the BLYP+vdW model, and from 2.88 Å to 2.83 Å for the PBE+vdW model. In both cases, the first minimum depth increases, resulting in the decrease of n TeTe . Another remarkable feature is the increase in the intensity of the second peak that shifts to low r values. The second peaks of g GeTe (r) and g TeTe (r) are associated to the correlations between Te atoms forming Ge-Te-Te and Te-Te-Te (or longer) chains. The dispersion forces lead to a reduction of the Te-Te inter-atomic distances. For this specific pair correlation function, the BLYP+vdW model is slightly more structured on extended range order distances, as shown by the moderate growth of a peak around 7.7Å in the g TeTe (r). Albeit limited, this peak contributes 131 4.3. Impact of the dispersion forces to improving the intensity of its homologue in the total g(r), bringing theory in better agreement with the experiment.

Focusing on the coordination numbers of Ge (n Ge ) and Te (n Te ) and on their changes when the vdW corrections are applied, they decrease from 4.53 to 4.34 and from 3.78 to 3.51 with PBE. With BLYP the changes are from 4.14 to 3.99 and from 2.56 to 2.31. It appears that the resulting coordination numbers are reduced and, as a consequence, the total coordination numbers gets closer to the one anticipated by the N-8 rule ( n = 2.7). This holds despite the fact that the total coordination number in the PBE case remains quite far from the N-8 rule value (3.68 against 2.7).

Looking at the topology of each system, we provide in table [4.2] the percentages of Ge atoms involved in homopolar bonds, edge-sharing and corner-sharing tetrahedra as well as the percentage of Te atoms involved in homopolar bonds. Concerning the Ge atoms, the subdivision proposed in table [4.2] is built on the existence of three classes (homopolar, ES and CS) and on the assumption that a given Ge atom cannot belong to more than one of them. The above numbers are consistent with the behavior of the partial pair correlation functions. Most remarkable is the PBE case, since the account of vdW forces largely increases (decreases) the number of corner sharing (CS) (edge sharing (ES)) connections. This is visible in the shape of the corresponding pair correlation function, in which the double peak pattern in between 3.5 Å and 4.5 Å is replaced by a bell shape centered around 3.8 Å.

N GeGe N TeTe N CS Ge N ES
A local environment analysis of the l-fold Ge and l-fold Te atoms are displayed in fig. [4.3]. The account of vdW forces leads to an enhanced structural organisation for both the XC functionals, through an increase of the fractions of 4-fold Ge and 2-fold Se. However, the network topologies are substantially different. In the PBE case, a sizeable fraction of Ge over-coordinated atoms do exist. This analysis will be further discussed when comparing the performances of the XC functional (next section).

Existing results for the Ge 15 Te 85 melt show that the account of the dispersion forces do not affect the angular distributions and the topology of the Ge 15 Te 85 melt [11]. In our case (Fig. 4.4), the effects of the dispersion forces are very small in the intra-polyhedra (Te-Ge-Te) bond angles distribution. A main peak at angles near the tetrahedral one (109˚) occurs for both models (BLYP+vdW and PBE+vdW). In contrast, for both PBE and BLYP, the shape of the inter-polyhedra bond angle distribution exhibits sizable changes. When turning on the vdW forces, the PBE main peak shifts from 82˚to 90˚. This shift is consistent with the increase of the fraction of CS configurations found for this model. The BLYP BAD(Te-Ge-Te) narrows and slightly increases in intensity. In contrast to the changes affecting the Te-Te inter-atomic distances and correlations, the BAD(Te-Te-Te) exhibits small changes. While the PBE+vdW BAD(Te-Te-Te) is less affected by the account of the dispersion forces, the BLYP+vdW main peak shifts slightly to higher values.

Impact of the XC functional

In this section, we compare the performances of BLYP and PBE exchange correlation fucntionals. In what follows, only models with account of vdW corrections will be considered.

Comparison between the simulated and the measured total structure factor s(k) and total pair correlation function g(r) is displayed in Fig. BLYP is superior in reproducing the position of the first peak in the total structure factor, although the first peak intensity is overstimated. The position of the second peak is accurately reproduced for both functionals. For k larger than 4 Å -1 , both functionals follow closely the experimental pattern, some (small) differences being noticeable. Turning to the total pair correlation function, BLYP gives better results since the intensities and positions of the first minimum and of the second peak located at 3.3 Å and 4.2 Å, respectively, are in better agreement with experiments. Fig. [4.6] provides the comparison between the partial g(r)s within the two models. g GeGe exhibits the most important differences between the two cases. When the dispersion forces are used, BLYP features a much higher intensity in the first peak. In addition to a deeper minimum, the same tendency is found in the case of g GeTe (r). Looking at g TeTe (r), the behavior of the first peak (less intense and narrow in the BLYP case) is compatible with a reduction in the number of homopolar Te-Te bonds. As displayed in table [6.1], BLYP with account of vdW corrections The choice of the XC scheme has a crucial effect on the amorphous structure. As shown in Fig. [4.7], the percentages of l-fold atoms found by using BLYP+vdW and PBE+vdW are substantially different. Both favor 4 coordinated Ge atoms as the predominant contribution. However, PBE+vdW contains 3-5-and 6-fold Ge atoms in a larger amount. For Se, the same occurs for 3-4-and 5-fold atoms. Figure [4.8] shows a comparison between the angular distributions with changing the XC framework. In the inter polyhedra BAD(Ge-Te-Ge), the main peak is found to be centered around ∼ 90˚. The BAD(Te-Ge-Te) changes drastically from PBE to BLYP, the main peak shifting to higher value (108 ˚) with the BLYP functional. This has to be correlated with the increase of the fraction of tetrahedral configurations. The small peak near 180˚of the PBE+vdW model is consistent with the existence of a small fraction of 5-fold Ge, as it will be explained when invoking the analysis of the local order parameter (later in this section). Indeed, peaks around 180˚are a sign of octahedral-like structures. The PBE+vdW (Te-Te-Te) BAD shows two mains peaks centered around 90˚and 180˚respectively, while in the BLYP+vdW model only one peak appear. These peaks are related to specific topologies due to the existence of Te chains (i.e. at least 3 Te atoms, one of them having 2 Te atoms as nearest neighbors) in both models. The 180˚angle is suggestive of a more important octahedral contribution (see what follows).

It is of interest to ascertain whether or not there is a coexistence of tetrahedral and octahedral coordinations in our models. To this end, we provide the distribution of the local order parameter (Fig. [4.9]). The occurrence of an intense peak near q=1 indicates the predominant tetrahedral arrangement in both models. The peak associated to the PBE+vdW model has a broad distribution including values in the range [0.8-1] that reflects the existence of GeTe 4 tetrahedra and the extended range of angles around the main peak in the (Te-Ge-Te)BAD. In the PBE+vdW model, the peaks in between q=0.2 and q=0.8 reflect the occurrence of different defective octahedral environments. In addition, the small peak occurring near q=0 firmly indicates the presence of a small fraction of octahedra. The BLYP+vdW model shows a more limited structural variety and only a broad peak in between q=0.2 and q=0.6. This peak is due to some defective octahedra, even though the peak position cannot be attributed to any specific arrangement.

In conclusion, the topology of glassy GeTe 4 strongly depends on the exchange-correlation functional adopted. The BLYP model is characterized by a tetrahedral arrangement with the presence of some defects. In the PBE model there is a coexistence of tetrahedral and octahedral arrangements.

Models annealed at T=300 K

We can also asses the effects of the exchange-correlation functional and vdW corrections by considering models studied at T=300 K without thermal cycle (i.e. without raising the temperature to achieve significant diffusion followed by a further quench T=300 K). After relaxing the 138 PBE model without vdW corrections at T=300 K, three independent models were generated, the three having in common the same initial conditions. In the first, we changed the functional from PBE to BLYP (model BLYP@300K). In the second, we changed the functional as above and we also activated the vdW corrections (BLYP+vdW@300K). In the third, we activated the vdW corrections within the PBE framework (PBE+vdW@300K). These three models were annealed at 300K for 50ps, 46ps, and 30ps respectively. The partial pair correlation functions are displayed in figure [4.10]. The account of vdW leads to a more structured shape of the pair correlation functions. This effect is not dependent on the length (or the occurrence) of a thermal cycle. Overall, the peak intensities are higher within BLYP. In addition to the shorter Ge -Ge bond length, the BLYP Chapter 4. Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4 140 model exhibits a deeper and better defined first minimum than PBE. The maximum at 6.5 Å for the PBE model shifts leftward by 1Å with BLYP. Looking at g TeTe (r), BLYP+vdW@300K provides a deeper first minimum and a higher second peak intensity than BLYP@300K and PBE+vdW@300K. Worth to notice is the flattening out of the peak around 6 Å found in BLYP+vdW@300K for g TeTe (r). 4.3 r i j : position of the first maximum of the pair distribution functions g i j ), partial coordination numbers n i j (integration of g i j up to the first minimum), coordination numbers n i and the mean coordination number n.

r GeGe (Å) r GeTe (Å)
The interatomic distances Ge -Ge and Te -Te found from the annealed models are different from those found by employing a thermal cycle. In contrast, the Ge -Te interatomic distances remain unchanged. As a consequence of the behavior of the Ge -Ge and Te -Te interatomic distances, higher coordination numbers are recorded. The BLYP@300K model and the BLYP+vdW@300K model give a mean coordination numbers smaller than the PBE@300K one.

An interesting observation concerns the effects of the vdW forces on the structure of amorphous GeTe 4 and the interplay with the characteristic times of relaxation of the glassy structure when changes in a given model are introduced. From the results of the present section, it appears that the structural changes occurring upon inclusion of the vdW correction manifest themselves on a very reduced time scale, i.e. without need of an extended thermal cycle. In other words, the impact of the vdW forces is readily observable on a relatively short trajectory at T=300 K. This is not true in the case of a change corresponding to the replacement of a given XC functional with another. For these situations, a thermal cycle is of paramount importance. To substantiate this assertion, we show (figure [4.11]) that the original and the new bond angle distributions obtained by changing the functional are very similar. Te -Ge (bottom panel). T=300K, Simulation BLYP@300K (black lines), simulation BLYP+vdW@300K (red lines) and simulation PBE+vdW@300K (green lines).

Chapter 4. Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4 142

Figure [4.12] displays the nearest neighbors analyses for the three models (BLYP@300K, BLYP@300K+vdW, PBE@300K+vdW). All networks feature a tetrahedral arrangements with a sizable fraction of 3-fold and 5-fold Ge atoms. The Te atoms reproduce the PBE+vdW@300K structural diversity with a mixture of 2-, 3-and 4-fold configurations. Due to the limited sampling of the configurations space, the BLYP results are quite similar to the PBE ones. However, one notices a small reduction of the 3-and 4-fold Te and 3-and 5-fold Ge. This trend is enhanced when activating the dispersion forces (BLYP+vdW@300K), confirming the BLYP tendency to promote the tetrahedral environment. 12 The fraction of l-fold Ge and Te atoms (l=1, 2, 3, 4, 5 or 6) as obtained from BLYP@300K (black color), BLYP+vdW@300K (red color) and PBE+vdW@300K (green color) models.

Part III Ternary Phase change materials

for storage able to overcome the intrinsic physical limitations and scaling problems when reducing devices dimensionality. Four major problems can be mentioned in the currently used memory devices. A memory device has to ensure a high data retention for a long period, this feature being related to the density of the used material. In addition, many of the used alloys suffer from phase segregation after some write/erase cycles which undermine the use of these device for long term storage. The third related problem is the slow switching speed from bit "1" to bit "0" and finally the high energy consumption required to write, read and erase data. In this context, the phase change materials (PCM) represent an excellent promising candidate to replace and fulfill the technological needs and built the commonly called "universal memory".

In the following we will sketch the physics related to the PCM in relation to their practical implementation. For an exhaustive reviews on the PCM technology, one can refer to [2,3]. PCM devices are built from chalcogenide alloys of the VI group of the periodic table. The concept of PCM devices is based on particular properties of chalcogenides that show reversible transition between two stable states of matter. The two states are characterised by resistance and/or optical contrast depending on whether the chalcogenides is amorphous or crystalline. The discovery of the phase changing electrical characteristics goes back to Alan Tower Waterman, during his studies of thermoionic emission of MoS 2 in 1917 [4]. Waterman studies showed that the MoS 2 exist in two forms, a high and low resistance states [5]. The progressive conductivity changes upon heating are a prominent feature of the phase change behavior in chalcogenide materials. Later, in 1968, the reversible electrical switching phenomena in disordered chalcogenides structures upon change in temperatures were investigated by S. R. Ovshinsky, considered as the father of PCM materials [6]. Ovshinsky found a practical implementation in the non-volatile PCM in 1970 [7,8]. Even though the PCM were discovered in the 70's, viable prototypes were made available only in the last years.

In the following, we will give hints on the reversible phase transition and we shall remind the major achievements on studying the ovonic systems.

Reversible phase transition

Phase change materials are chalcogenides that exist in amorphous and crystalline stable phases. The remarkable differences on the electronic and optical properties of both phases are due to the change in the bonding upon the crystallization of the amorphous phase. The change in bonding and the related mechanisms induce the change of the physical properties of these materials. 148

Reversible phase transition

To be considered as a good PCM, a material should ensure the optimal compromise between Figure 5.1 Phase transition from amorphous to crystalline state in a PCM material. A long laser pulse (set operation) with moderate-intensity beam heats the PCM spot above its glass-transition-temperature leading to a high mobility regime that converts to the crystalline energetically favorable phase. A short laser pulse (reset operation) with high-intensity beam heats the crystalline spot above its melting point. By cooling via the contact with the surrounding amorphous, this crystalline spot converts to the amorphous state. The crystalline state is an ordered state characterized by a low resistance compared to the amorphous state, thus a measure of the resistance allows to read the stored information, whether binary "0" or "1". Figure modified from the original figure 1 in [9].

the following properties:

• High speed phase transition.

• High stability of amorphous and crystalline phase for a long period of time at the operating environmental conditions.

• Large contrast in optical properties and/or resistance change.

• Large cycle number of reversible transitions.

• High chemical stability.

• High laser light absorption.

Figure [5.1] illustrates the principle operations of a memory devices based on PCM. For further details and explanation we refer to refs. [9][10][11] and included references. Given that the system is in the amorphous state at ambient temperature (considered as the "reset" state equivalent to '0' binary), we heat locally the PCM using a low intensity laser beam (voltage) to the crystallization temperature (T x ) (above glass-transition temperature t g ). This allows the atoms to reach a high mobility without melting and revert for enough long heating time to the energetically favorable state. The crystalline spot (considered as the "set" state equivalent to '1' binary) is an ordered state and it is characterized by a low electrical resistivity or/and a high optical reflectivity. As such, it could be distinguished from the surrounding amorphous alloy. In this way we write a "1" bit. To erase the stored information, one should heat the crystalline spot using a high intensity laser beam above the melting point (T m ), leading to the formation of a liquid state. By contact with the surrounding lower temperature amorphous regions, the liquid will be rapidly cooled/quenched to the amorphous state. Reading the stored information is simply a measure of the resistance or the optical reflectivity over the alloy.

Currently studied PCM

The first practical demonstration of the switching phenomenon was carried out by Ovshinsky [6]. The material used was based on the Te-based eutectic alloys, doped with S, Ge and As. Although these materials showed the possibility of a real application on a memory devices, they suffered from a very slow crystallization speed (of the order of microseconds) that limits their use. A faster crystallisation was found for GeTe alloys [12] and Ge 11 Te 60 Sn 4 Au 25 [13], the study of these materials triggered the discovery of good PCM along the GeTe -Se 2 Te 3 tie line. Currently, Ge 2 Sb 2 Te 5 [14][15][16] is considered the benchmark of PCM and it has been used in many commercial devices. A second generation of PCM are those based on Sb 2 Te alloys doped by elements such as In, Ag, Ge [15][16][17]. Recently a third generation of materials of Ge-doped Sb material has been proposed [18,19], and chapter 10 in [2].

Table 5.1 Chronology of phase-change materials developments [2,20].

The secret behind PCM materials

Why only chalcogenides could be used to built optical memories? More precisely, how can PCM's show a high crystallization speed and at the same time a good thermal stability of the amorphous phase? On the one hand, the underlying reason is related to the electronic structure of chalcogenides. Te and Se can make two fold coordinated bonds with nearby atoms, usually weakly bonded by van der Waals forces (e.g. case Te-Te-Te). This kind of bonds are easily broken when heated to the melt state, since the melting temperature is usually low. On the other hand, one common feature of PCM is the high symmetry in the atomic distribution. For example, if we take the extensively studied tie line GeTe -Sb 2 Te 3 , the majority of the compounds belonging to this region do exist in a rocksalt structure. For instance, figure [5.2] shows the structure of Ge 2 Sb 2 Te 5 . It can be expressed as simple cubic blocks were the Te occupy one rocksalt sublattice and the rest are randomly filled by Ge and Sb, with the coexistence of large distortions and voids even at low temperature. Also superior PCM like Ag 3.4 In 3.7 Sb 76.4 Te 16.5 shows high symmetry, and interestingly with a random distribution of all atoms. For all these materials, the molten high temperature regime is characterized by a very similar distorted network (cubic, rocksalt) and only small atomic rearrangements leads to an ordered crystalline phase. This feature is related to a particular thermal property, since many of PCMs show a very and amorphous (right) states. Ref [24].

Using reverse Monte Carlo simulation, Kohara and Co [25] simulated the crystallization of amorphous Ge 2 Sb 2 Te 5 and GeTe. This study showed that the main features of the crystalline structure such as the even membered rings and bond-angles centered around 90˚remain in the amorphous state with small distortions. Thus, the crystalline structure is easily and rapidly recovered upon small atomic rearrangements with heating (see Fig. [5.4]).

Figure 5.4 Schematic presentation of the possible crystallization process in amorphous

Ge 2 Sb 2 Te 5 and GeTe. Ring structure transformation in the phase changes crystal-liquid-amorphous. Ref [25].

Hegedus and Elliott showed that the fast crystallisation of Ge 2 Sb 2 Te 5 is associated to the formation of a chemically ordered medium-range ordered (planar four-fold rings, Figure [5.5]) structures driven by the intrinsic electronic structure of the atoms and their bonding properties, which facilitates nucleation and growth of the crystal seeds [28]. 5.5 Ga -Sb -Te as a good PCM??

Recently, interest has been devoted to discover new phase-change materials with outstanding performances. In particular, the research was targeted to enhance the properties of the commonly used Ge 2 Sb 2 Te 5 that showed a proven adequacy for NOR flash memory applications. This material showed a growth dominated crystallization mechanism, while suffering from low thermal stability drawbacks. Nowadays, applications such as automotive require a material with short crystallization time, high crystallization temperature and large resistivity contrast between the amorphous and crystalline states.

In 2005, C. M. Lee introduced a new chalcogenide materials based on Ga-Sb-Te films [33,34]. It is known that the Ga-Sb phase change material and doped composition exhibits high crystallization properties very similar to doped Sb-Te materials [35], but in contrast to the latter Ga-Sb shows an enhanced amorphous phase stability due to the presence of Ga [36]. In his study, C.M. Lee showed that recording disks based on GaSbTe media have demonstrated a superior writing/erasing performances. In addition, the crystallization rate has been highly enhanced when doping GaSb into Sb 7 Te 3 or Sb 8 Te 3 tie-lines [33,34]. The characterisation of many compositions along these tie lines showed that the Sb-rich compositions feature a lower melting temperature and a high T x /T m ratio, good amorphization and a good thermal stability. Many compositions showed good characteristics to be a potential PCM, Ga 2 Sb 5 Te 3 [33,37], Ga 26.4 Sb 65.2 Te 8.4 [38,39] and Ga 3 Sb 8 Te [38,40].

The crystallization behavior of binary Ga-Sb materials with different concentrations of Ga were studied by Cheng and Co. [41,42]. The stoichiometric GaSb has shown a very fast 154 155 5.5. Ga -Sb -Te as a good PCM?? crystallization speed, with a low difference in the contrast between the amorphous and the crystalline phase. The off-stoichiometric alloys feature a high phase stability (crystallization temperature higher than Ge 2 Sb 2 Te 5 (T x )). An interesting point raised in this study is that doping the GaSb alloy by Sb and Te leads to increases T x . Figure [5.6] show the explored compositions in the Ga-Sb-Te phase phase diagram in the literature. Among all the ternary Ga-Sb-Te compositions, the Ga 4 Sb 6 Te 3 [42] composition is partic- Crystallization temperatures T x as a function of composition in the Ga-Sb-Te phase phase diagram. Ref [42] ularly worth mentioning. This melt quenched alloy showed a fast crystallization time (190 ns) and a high crystallization temperature 271˚C in addition to a high electrical contrast R a /R c = 2.7 × 10 4 and low melting temperature (∼ 600k). In one hand, this characteristics are by far better than the commonly used GST, which has t x = 150 ˚C and in the other hand make comparable to the proposed Ga 2 Sb 5 Te 3 [37] as promising candidate for phase-change memory.

A very recent paper by Kao & Co [43] demonstrated that another composition in the GaSbTe phase diagram, the Ga 18 Te 12 Sb 70 , is a very good PCM candidate. It shows a high T x = 248 ˚C and a temperature corresponding to 10-year-failure, 183 ˚C. In addition, a reset and set speeds of 20 ns were recorded with two-order of magnitude change in electrical resistance. In general, this work concluded that the composition ranges Ga 18 ∼ 25 at. %, Te 8 ∼ 12 at. % and Sb 67 ∼ 70 at. % are suitable to ensure thermal stability with T x > 240 ˚C and an activation energy of crystallization > 5.0 eV. Despite the promising features of GaSbTe alloys that makes them an excellent candidate for universal phase-change memory, the understanding of their atomic scale structures remain very poor. No first-principle studies on GaSbTe alloy were available when we started working on this topic. Our results are intended to fill this gap, by providing the structure and the electronic properties of one of the most intruiging PCM alloys, Ga 4 Sb 6 Te 3 .

STRUCTURE OF THE AMORPHOUS

Ga 4 Sb 6 Te 3 VIA FIRST-PRINCIPLES MOLECULAR DYNAMICS

Let's Start

This work was carried out in collaboration with the team of Professor Marco Bernasconi in Milan, University of Milano Bicocca. This collaboration led us to exploit and compare the standard Car-Parrinello MD approach and a novel approach derived from this one, termed the second generation of this method (SGCPMD). Very shortly, one can say that the main difference between the two methods lies in the non self-consistent optimization of the wave functions during the dynamics. This method was proposed by Kuhne & Co [1] and implemented in the CP2k suite of programs [2,3]. In contrast to CPMD, SGCPMD allows the use of large integration time steps which lead to a slightly dissipative dynamics of the type -γ D ṘI , where R I is the ionic coordinates and γ D is an intrinsic friction coefficient. Further details could be found in [1]. In order to study the amorphous Ga 4 Sb 6 Te 3 , we adopted a generalized gradient approximation (GGA) for the exchange and correlation part of the total energy. In particular, the BLYP scheme was used since a recent study [4] showed that BLYP leads to better agreement with experimental structural properties than PBE in the case of liquid and amorphous Ga 2 Sb 2 Te 5 . Martins-Troullier norm conserving pseudopotentials are used for the core-valence interactions. For CPMD calculations, the wave functions were expanded in plane waves basis set at the Γ point of the supercell. The calculations have been performed with the non linear core correc-tions (NLCC) for Ga and a plane waves cutoff of 60 Ry [5,6]. An integration time step of 0.12 fs and a electronic fictitious mass of 600 (a.u) were used. In order to implement the Quickstep scheme [7] (see appendix D for details) (which is a prerequisite to the application of the second generation Car-Parrinello dynamics, at least for all implementations carried out by the team of Prof. Bernasconi), the Kohn-Sham orbitals were expanded in a triple-zeta-valence plus polarization (TZVP) Gaussian-type basis set and the charge density was expanded in a plane wave basis set with a cutoff of 100 Ry to efficiently solve the Poisson equation within the Quickstep scheme [7]. Brillouin zone integration was restricted to the supercell Γ points. A time step of 2 fs was used for the sample preparation and data collection simulations. A temperature control is implemented for both ionic and electronic degrees of freedom by using the Nosé-Hoover thermostats. In order to generate an amorphous model of Ga 4 Sb 6 Te 3 , two different initial configurations were considered: First, a small model of 117 atoms was generated starting from an amorphous GeSe 2 system by changing the atomic identity of the atoms. Ge atoms were replaced by Ga atoms, and Se atoms were randomly replaced by Sb and Te atoms. The size effect on the structure and the local order in the liquid and amorphous Ga 2 Sb 2 Te 5 was shown to be not significant [4]. Since a system as small as 117 atoms is not necessarily expected to provide a detailed description of the crystallization mechanism, larger systems are needed. Known that the Ga 4 Sb 6 Te 3 crystallize in a rhombohedral phase with lattice constant a=b=4.28 Å and c=17.04 Å [8]. A fairly large network was generated with a rhombohedral arrangement using VESTA code in which the atomic plans were randomly distributed. From this model we produced a cube of 299 atoms. Figure [6.1] shows some snapshots of the initial configuration.

Fixing the liquid density

Since the experimental density is unknown, the initial density of the two models was adjusted to 0.337391 Å -3 (6 g/cm 3 ), intermediate density between binary GaSb, [9] and ternary Ge 2 Sb 2 Te 5 [10]. In order to find and calibrate the system density, the N=117 model was simulated using CPMD. Initially the system was heated up to T=2000 K in 200 ps, equilibrated and further annealed for ∼ 25 ps, quenched to the liquid state at T=900 K and equilibrated for 5ps. By computing the stress tensor on the final configuration, we calculated ∼ 2.77 GPa for the internal stress. In order to release this stress, the system density was relaxed by increasing the cell side lengths and computing the stress tensor. The density that minimises the internal stress was found to be 0.030755 Å -3 (alat=29.5 a.u). In addition to this, a run on the isobaric isothermic ensemble (NPT) was performed for more than 20 ps at T=900 K (without NLCC for Ga atom). Figure [6.2] shows the evolution of the instantaneous internal stress and the lattice parameter during the run. The cell parameter was equilibrated around 29.47 atomic units and the internal stress fluctuates around 0.37 GPa, consistently with the "hand made" calibration of the cell. The density of the liquid state was fixed to the value obtained from the NPT run, 0.031 Å -3 (alat = 29.5 a.u ) at T=900 K and the system was further annealed for over 25ps in the NVT ensemble. Figure [6.3] shows the evolution of the diffusion coefficients of Ga, Sb and Te atoms as calculated from the mean-square displacement of a given species at time t. The inset shows the time dependence of the mean square displacement of different species. The diffusion coefficients of Ga, Sb and Te are respectively: 0.24 10 -5 cm 2 /s, 0.16 10 -5 cm 2 /s and 0.14 10 -5 cm 2 /s.

Liquid state and shortly quenched amorphous

The liquid state has been quenched from T= 900 K down to T=300 K on 100 ps. At T=300 K, we relaxed and calibrated the density to 0.0278 Å -3 for which the internal stress is negligible. Data were collected for 25 ps at T= 300 K. Figure [6.4] order typical of a liquid. When quenched, all the peaks sharpen and increase in intensity. The amorphous state shows an extended order in the region [5][6][7][8][9][10] Å where many new peaks grow up with decreasing the temperature. The first peak in g GaGa (r) moves from 2.70 Å at T=900 K to 2.55 Å at T=300 K. The first minimum position moves toward low r values, leading to a reduction of the partial coordination number n GaGa and to a reduction of the fraction of the homopolar bonds (see table [6.1]). Compared to the behavior at T=900 K, g GaSb (r) and g GaTe (r) feature a higher first peak at 2.77 Å at T= 300 K. This is associated to the increase of the chemical order with the quench. The first minimum is well defined at ambient temperature, this fact is reflected by the decrease of the associated partial coordination numbers. In addition, two peaks arise at 4.25 Å and 6.69 Å corresponding to the appearance of a second and third neighbor shells.

Amorphous state with long quench from CPMD and SGCPMD

It is worth noticing that the Ga -Ga bond lengths found in this model are larger than those found in crystalline GaTe (r Ga-Ga = 2.43 Å ) [11], crystalline orthorombic Ga [12] (r Ga-Ga = 2.48 Å ) and in Ga/Sb films [13] (r Ga-Ga = 2.47 Å ). g SbSb (r) and g SbTe (r) are characterized by an intense first peak at ∼ 2.9 Å, reflecting the abundance of the Sb -Sb and Sb -Te bonds. Focusing on the environment of the Sb atoms, in the region comprised between 3.5 Å and 10 Å, there are three peaks indicating the extended range order and the preference to be linked to Sb as first neighbor and Sb or Te as second neighbor. At T=900 K, the pair correlation function g TeTe (r) has a high broad first peak around 4.13 Å. This peak moves to 4.2 Å upon quenching. While homopolar bonds are not present in crystalline GaTe [11], a few homopolar bonds are found in the liquid state (cutoff distance of 3.2 Å) and disappear with quenching. The g TeTe (r) distribution has large statistical fluctuations between 6 Å and 9 Å, mostly ascribed to the small number of Te atoms. All the atomic coordination numbers reduce with the quench. In striking contrast with amorphous GaSb [13], where Ga is fourfold, the Ga atoms appear to be mostly threefold coordinated. The Sb coordination number compares well to the coordination number of Sb found in the case of amorphous GaSb 7 [13]. However, this number is smaller than the one calculated for amorphous GaSb [13]. This disagreement is somewhat intriguing and calls for further calculations in the attempt to elucidate its nature. For this reason, we have undergone an additional set of calculations by using a longer quenching schedule.

Amorphous state with long quench from CPMD and SGCPMD

By terming model A the one described and employed in the previous section, we now introduce a larger system made of N=299 atoms and longer quenching schedule. Accordingly, we can define three models as follows:

• Model A-117 atoms, quenched from the melt state during 100 ps with CPMD.

• Model B-117 atoms, quenched from the melt state during 300 ps with CPMD.

• Model C-299 atoms, quenched from the melt state during 300 ps with CP2K.

Long quenched models B and C were generated by using the plane waves framework (code CPMD, Model B) and the combined plane waves/gaussian scheme ( As noticeable in table [6.2], the Ga-Ga interatomic distances calculated within model B and model C are different. In order to assess whether this can be ascribed to the nature of the pseudopotentials, we considered a supercell of crystalline GaTe (N = 240) [11] by using two different families of pseudopotentials (GTH and MT). The final results are as follows:

distance (Å) CPMD-MT CPMD-GTH exp. As a further test, we relaxed the structure of the crystal by using the package Quantum Espresso (8x8x8 k point mesh, 60 Ry cutoff) that is based on plane waves and pseudopotentials much in the same spirit as CPMD. In addition, this code is more suited to study crystals since the k-point sampling can be easily exploited. By taking the GTH pseudos (please note that MT pseudo are not implemented in this code) we obtained the same interatomic distances as the ones provided by the CPMD with GTH pseudos. The experimental Ga -Ga distance for crystalline GaTe lies in between the one found from CPMD-GTH and CPMD-MT. Differences 6.2. Amorphous state with long quench from CPMD and SGCPMD could be ascribed to the nature of the used pseudopotential, therefore we conclude that CPMDbased distances can be taken as reliable.

The first and second peaks in g SbSb (r) and g SbTe (r) obtained with model B and model C are in very good agreement. Taken with respect to Model A, a little decrease in the intensity of the Sb-Te first peak is reflected by the increase in the corresponding partial coordination number. In all of the models considered (Model A, Model B and Model C), no homopolar Te-Te bonds were found. Comparing the results for the two distinct quench rates, it turns out that only a little effect can be observed in the position and the intensity of the main peak of g TeTe (r). The resulting coordination numbers are given in table [6.4] The coordination numbers are strongly affected by changes in the quenching schedule. Overall, they appear to increase for the longer quench. For instance, the Ga coordination increases to reach a value close to 4 for the C model. The values for Sb and Te agree well for model B and model C. As a consequence, the mean coordination numbers increase as a result of quenching. These results confirm that 300 ps are needed to reach a coordination number very close to 4 for Ga. This is consistent with the value reported for the crystalline structure of GaTe (n Ga = 4) [11,14] and the near eutectic composition amorphous GaSb 7 (n Ga = 3.98) [13]. Moreover, our results are in agreement with those reported in Ref. [13], for which Sb was found almost three fold (n Sb = 3.19). Te atoms are bonded in average to 2.8 atoms, in agreement with the crystalline GaTe structure.

For our quenched amorphous models, the bond angle distributions (BAD) are displayed in figure [6.7].

The main peak for BAD (Sb-Ga-Sb) is close to the tetrahedral value 109.5˚. The small bump existing for model A near 91˚is seen to disappear after quenching. The main peaks of BAD (Te-Ga-Te), BAD (Ga-Sb-Ga) and BAD (Ga-Te-Ga) are centered near 100˚, that is inbetween the tetrahedral and the octahedral values. The Sb centered configuration shows a broad distribution of angles. In particular, for Te-Sb-Te, a first peak is located around 95˚for model B and C

The Ga environment shows an equal richness in Te and Sb in agreement with the partial coordination numbers (n GaSb = 1.78 and n GaTe = 1.80). The Ga environment consists of Ga -Sb n Te 4-n (Model B: 64.2% and Model C: 59.09%, n=0, 1, 2, 3 or 4) and Ga -(Ga -Sb n Te 3-n ) (Model B: 18.10% and Model C: 25.88%, n =0, 1, 2 or 3). This arrangement is in agreement with the main peak found in the Te-Ga-Te and Sb-Ga-Sb bond angle distributions and the local order parameter (see appendix C for a definition) (see Fig. [6.8]).

The local order parameter has a main broad peak centered around 0.95 indicating a predominant and yet not perfect tetrahedral arrangement. Due to the very rapid quench adopted to produce Model A, the fraction of fourfold coordinated Ga atoms is limited to ∼65 %. This tetrahedral geometry around Ga is remindful of that found on GaSb alloys [13], but in our case many of the Sb atoms are replaced by Te atoms. Ga-Ga-XX motifs exist not only GaSb alloys but also in the GaTe crystal [11,14]. This latter is made by Ga -GaTe 3 groups tetrahedrally arranged.

By comparing the coordination numbers relative to Sb, one realizes that n SbSb = 1.92 is larger than n SbGa = 1.19, moreover the fraction of Te around Sb atoms is even smaller (n SbTe = 0.19), this being a possible sign of phase segregation. When looking at the geometrical arrangement around Sb, two main arrangements could be distinguished: a dominating trigonal Sb -Ga n Sb 3-n (Model B: 50.91%and Model C: 55.88%), n =0, 1, 2 or 3 and a tetrahedral one Sb -Ga n Sb 4-n (Model B: 27.53%and Model C: 27.04%, n =0, 1, 2, 3 or 4). The trigonal Sb centered motifs are the main features detected in Sb-rich GaSb 7 film and coexist with a small fraction (∼17%) of 4-fold Sb [13]. This features have been also found in Ge 2 Sb 2 Te 5 [15]. The local order parameter for all the species are displayed in figure [6.8]. The Sb local order parameter exhibits a very broad distribution due to the variety of Sb fourfold environments. Te is mostly threefold coordinated, with a predominant motif Te -Ga n Sb 3-n (Model B: 70.61% and Model C: 68.40%, n =0, 1, 2 or 3). The above motif is also found in the GaTe crystal together with Te -Ga 3 blocks. In the Te coordination environment one also notices the presence of Te -Ga n Sb 2-n planar features (Model B: 27.73% and Model C: 29.27%, n =0, 1 or 2).

Overall, the structure is a mixture of pyramidal, tetrahedral and planar motifs which can not be described by the local order parameter as defined for the tetrahedral and octahedral geometries. A new definition accounting for motifs ranging from from planar (or pyramidal) to tetrahedral structure is required. For Ga two motifs coexist: the dominant tetrahedral geome- tries and some pyramidal like structures. Inversely, the Sb environment is mainly trigonal with the existence of tetrahedral geometries. Finally, Te is found in a defective trigonal (pyramidal) structure, with some plan structures with obtuse angles. Definition and computational details to get the void distributions are provided in appendix C. In PCMs the voids play a crucial role in rapid crystallisation. Akola and Jones [16] proposed that the presence of voids on GST and GeTe favors the reorientation of four membered rings, leading to accelerating the crystallisation speed. Cavities are found in our GaST systems are larger than those recorded for GST and GeTe [16]. The presence of nanocavities might be the reason behind the high crystallisation speed in amorphous Ga 4 Sb 6 Te 3 . The occurence of such voids could be due to the presence of large size Sb clusters (see Fig. [6.20]), instead of being homogeneously distributed the Sb atoms form "compact" clusters. The larger cluster found in model C, being made from 97 Sb atoms. The topology described above can be considered as a possible sign of a phase separation. This kind of phenomenon is an important issue for applications in the area of phase change materials. [17][18][19]. Such a phase separation involving Sb and GaSb was claimed to be unavoidable during the crystallization in Ga 16 Sb 84 [20]. In our case, the Ga 4 Sb 6 Te 3 do not show such a phase separation but one can conclude that our system can be envisioned as (GaTe) 3 -GaSb 6 , in which Sb-Te bonds are rare and nano Sb clusters can be found (see Fig. [6.12]).

6.3. Models annealed at 300 K Overall, and for all partials, the three models yield very close results. All the first minima are now well discernable, and the positions of all peaks are reproduced within a small error. The first peak position (FPP), the first minimum position (FMP), and the partial coordination numbers, as computed from the PCFs, are displayed in the following When comparing these results to those of the previous sections, one can notice that when changing from MT pseudos (code CPMD) to the GTH pseudos used in the GPW/SGCPMD approach (in short referred to as CP2k), the interatomic Ga-Ga and Ga-Sb distances slightly shorten. The same distance increases in the opposite case, that is to say when a CPMD like relaxation begins from a configuration created in the CP2k framework. This behavior could be ascribed to the nature of the pseudopotentials and/or to the pseudpotential cutoffs employed. This issue deserves further investigations.

In contrast to the melt quenched models, all annealed models have a similar Ga -Ga distances, however the FMP associated to the g D GaGa remain shifted to higher r values. The lower intensity of this first peak is such that, by compensation, the partial coordination number from model D agree with model F. Model E and F are less affected by the annealing and the associated coordination numbers remain very similar to the original models (B and C respectively). Model D exhibits dramatic changes in the coordination numbers with respect to Model A, indeed the obtained structure is quite similar to the one of models E and model F. This result confirms that we can rely on the structures of models quenched during 300 ps as a legitimate description of this glass. In addition, Model D (stemming from the original Model A) gives a structural description quite similar to that of Model B and C. That is to say, the annealing at T=300 K has the effect of rendering the structure obtained for the small sample very close to those obtained for the larger samples.

General conclusions and openings

Binary and ternary chalcogenides are at the heart of most technological applications that we are using. Modelling their structures at the atomic scale level is a way to link the macroscopic physico-chemical properties to the microscopic features in the bonding. This work takes advantage of the available high performence supercomputers and the most developed atomic scale methods to provide a robust understanding of chalcogenide materials. The so called Car-Parrinello molecular dynamics was widely used to model amorphous structures. We resorted to the first generation of this technique ("standard CPMD") to model the GeSe 2 , GeSe 4 and GeTe 4 systems, while we combined CPMD to its second generation (SGCPMD) scheme to produce the Ga 4 Sb 6 Te 3 glassy structure. Our research led to a better understanding of the behavior of binary Ge x Se 1-x systems under pressure. Furthermore, in addition to the effects of vdW forces and the performances of exchange-correlation in modelling Te-base materials, we generated an amorphous structure of Ga 4 Sb 6 Te 3 for which we obtained an accurate description of the glassy network. In what follows we recall the main results of the present thesis:

In the case of GeSe 2 and GeSe 4 amorphous systems we obtained the following results:

☞ First, we worked on enhancing the ambient pressure model of glassy GeSe 2 . Investigation of the simulation methodology revealed that a residual pressure of ∼1 GPa alter the network topology and leads to discrepancies with results from neutron-diffraction experiments.

By releasing this pressure we were able to generate a model that reproduces to a large extent the experimentally measured topological features. For instance, the fractions of Ge atoms belonging to edge-sharing and corner-sharing configurations and the fractions of homopolar Ge and Se atoms. We also demonstrated that the account of residual pressure effects on the room temperature, equilibrium structure of the glassy GeSe 2 leads to a better agreement between atomic-scale models and experiments for this prototypical network-forming disordered material.

☞ Second, we generated independent glassy GeSe 2 systems under pressure in both NPT and NVT ensembles. Our equation of state reproduces the X-ray measured one and covers a range of pressure up to P∼16 GPa. The computed reciprocal space quantities, say the total and partial structure factors, showed a good agreement with isotopic substitution neutron diffraction experiments. The behavior of this glass under pressure was found to be substantially different from models described in the literature. The present work proved the existence of a structural phase transition from a chemically-ordered tetrahedral network to a pseudo-cubic network with an average coordination number n ∼3.9. This transition is described as a two stage transition. In the first stage, as the density is increased up to the threshold pressure of ∼8.5 GPa, the chemical order is maintained and the network topology is less affected. The mean coordination number remains constant and small fluctuations of the numbers of CS and ES tetrahedral units occur. In the second stage, as the density is increased beyond P∼8.5 GPa, the mean bond length and the mean coordination number increase steadily. The CS to ES ratio reduce monotonously and 4fold coordinated Ge atoms are found in combination with 5-and 6-fold coordinated Ge atoms. We demonstrated that homopolar bonds drive the structural transition. Finally, the network arranges in a pseudo cubic fashion where all the A-B-A bond angle distributions are centered around ∼90˚and ∼180˚. The transition mechanism is in contrast with others pertaining to the crystalline phase, where transitions are observed from an ambient pressure 2-D structure in which layers are formed from equal numbers of ES and CS tetrahedra, to 3-D structures of densely-packed CS tetrahedra. Also, the description of the transition mechanism proposed in this work is in contrast with the one proposed by Drabold and Mei where the densification mechanism was interpreted in term of continuous increase of the mean coordination number.

☞ We studied the GeSe 4 amorphous system using the same theoretical recipe as in the case of GeSe 2 . We demonstrated that the GeSe 4 network does not undergoe any structural transition as the pressure is increased. Indeed, the chemical order and the network topology remain very less affected for pressures up to P∼13.2 GPa. The mean coordination number as well as the different bond lengths remain constant. Thus, the densification does not alter the short range order. However, through the analysis of the real space properties (the partial pair correlation functions), we showed that the densification of this glass occurs through a rearrangement in the second shell of coordination and the higher order shells around Se. We elucidated the mechanisms behind the persistence of the chemical order under load. It lies on the Se chains, which acts as springs linking the GeSe 4 tetrahedra. When subject to pressure, Se chains bind and twist by reducing their internal Se -Se -Se angles by filling the inter-tetrahedral space by avoiding overcoordination of the atoms.

General conclusions and openings

In the case of GeTe 4 system we obtained the following results:

☞ This study addressed two main issues. The first issue is related to the effects of the vdW corrections on the simulated GeTe 4 glassy structure. We simulated models with and without account of the dispersion forces and with two different XC functionals. We demonstrated that, regardless the XC recipe used, vdW correlations enhance the short and intermediate range structural properties when turned on. The mean coordination number decreases and the tetrahedral arrangement is enhanced. A very small effect on the bond angular distribution was found. Overall, the simulated structure with account of vdW empirical corrections leads to better agreement with experiments.

☞ The second issue has to do with the understanding of the performances of XC functionals when it comes to describing the GeTe 4 glassy structure. We focused on the performances of BLYP and PBE exchange correlation functionals. When comparing the measured total structure factor and total pair correlation functions, it appeared that BLYP gives higher level of structural agreement. One should note that in the case of Ge -Ge partial pair correlation function, BLYP and PBE show large differences. These discrepancies reduces slightly when vdW is activated with both functionals. The topology of the GeTe 4 network as described by the two XC recipes shows substantial differences. BLYP favors a tetrahedral network (4-fold Ge and 2-fold Se) and only a small fractions of over coordinated atoms was found. In contrast, PBE seems to favor the coexistence of tetrahedral and octahedral environments. The PBE model contains a larger fractions of 5-and 6-fold Ge and 3-, 4-and 5-fold Se atoms. Consequently, the angular distributions show different shapes with peaks centered near 109˚for BLYP models and ∼ 80t o ∼ 100˚for PBE models. We confirmed the coexistence of the two mixed tetrahedral and octahedral environments through the calculation of the local order parameter distribution for PBE model.

☞ On the general point of benchmarking and assessing the performances of existing XC functionals for the treatment of Te-based systems, BLYP with account of the vdW dispersion forces provides the highest level of topological accuracy when compared to experimental results.

In the case of Ga 4 Sb 6 Te 3 system we obtained the following results:

This study was carried out in two steps:

☞ First, we generated a small model (N=117 atoms) of Ga 4 Sb 6 Te 3 by melt-quenching during 100 ps using the CPMD approach. The liquid and amorphous densities were fixed by relaxing the structures to reach near zero internal stress. By analysing the obtained amorphous structure we showed the existence of unusual topological features that could be related to the short quenching schedule. To clarify this point, a second model was melt-quenched to T=300 K during 300 ps. We demonstrated that the quenching schedule impacts the final network topology. The atomic coordination numbers of the 300 ps quenched amorphous Ga 4 Sb 6 Te 3 compares well with the one of amorphous GaSb and GaSb 7 alloys and the crystalline GaTe structure. In addition, less structural variety was found, confirming that at least 300 ps of quenching are needed to get reliable amorphous structre of Ga 4 Sb 6 Te 3 .

☞ In the second step we generated a larger amorphous model (N=229 atoms) by applying the SGCPMD approach and longer quenching schedule (300ps). Our results suggest the absence of size effects affecting the final structure. Indeed, a self consistent picture is obtained from both simulations (CPMD with N=117 and SGCPMD with N=299). The Ga atoms were found in a coordination close to four, while Sb presents two different types of local environments, three-fold coordinated octahedral-like and tetrahedral. Te can take a 2-fold and 3-fold coordinations. The structure is also similar to the bonding network of the crystalline phases of GaTe and GaSb, being composed prevalently by Ga atoms in a tetrahedral environment.

☞ We demonstrated the existence of very few Sb-Te bonds in the amorphous phase. Thus, our Ga 4 Sb 6 Te 3 amorphous is a combination of two binary systems GaTe and GaSb (say Ga 4 Sb 6 Te 3 ≡(GaTe) 3 -GaSb 6 ). We showed that the high concentration of Sb leads to a partial segregation of antimony forming clusters with structures that recall the metallic phase of crystalline Sb.

189

General conclusions and openings

Future goals

❅ For what concerns the binary Ge x Se 1-x systems, we are willing to take advantages of the large data base of systems that we produced under pressure. In particular, additional interesting physical properties could be explored. Precisely we can think of Raman, the electronic, and the infra-red properties.

❅ The present work provided an interesting description of the structure of one phase change material, the Ga 4 Sb 6 Te 3 . The electronic and vibrational properties of this PCM should be investigated to explain the nature of some specific bonding (e.g. Sb clusters). Since the rapid crystallisation is one of the crucial properties of PCMs, we are aiming at exploring the event of a spontaneous crystallization of the glassy state by annealing the amorphous structure at the crystallization temperature. The crystallization kinetics and mechanisms deserve also to be elucidated.

❅ On a more general level, the knowledge developed during this PhD project will allow to explore and investigate new materials of high technological pertinence. In particular, we are aiming at studying the substituent effect on the glass structure and physical properties of a new generation of far-infrared transmitting telluride materials. Among the proposed systems the GeSe 4 -GeTe 4 pseudo-binary composition line of the Ge-Se-Te ternary diagram is of particular interest and also Ge 20 Te 70 Se 10 . In addition, we are willing to extend our studies to compounds of interest for infrared photonics such as the germanium-based telluride glasses (Ge-Ga-Te and the Ge-Te-I ternary systems).

Appendices

A-General methodology to generate binary amorphous models under pressure

In this section, we shall remind the general practical implementation to generate our amorphous models. One of the most popular methods to generate amorphous systems is the melt-quench method. As a starting choice, one can fix the system density to the amorphous density. Through randomization at high temperature beyond the melting point, this favors atomic mobility and efficient sampling of the configurational space. Once the atomic mobility reaches levels typical of liquid diffusion (10 -5 cm 2 /s), one can cool down the system to the melting temperature T m . This regime is characterized by a linear variation of the mean square displacement with respect to the simulation time. To get the amorphous, one should quench the melt to the ambient temperature. In principle, this method is not devoid of drawbacks, since the quenching schedule could lead to different amorphous structures. In this work, we did not specifically address this issue, since it was demonstrated (in other research efforts of our team) that an extended relaxation at room temperature is more significant (and can lead to much better structures) than variations of the quench rate within 1 order of magnitude. In principle, as shown in figure [6.15], one could revert back to a crystal from the liquid if the quench rate is low enough. However, this possibility is quite remote for most disordered system, because of intrinsic time limitations of our simulations. To summarize, we used the lowest quench rates compatible with the expected affordable computational effort and we focused on the structural relaxation of the disordered systems close to room temperature. This structural relaxation is more important than any expected structural difference due to different quench rates. In most cases, we selected as a starting configuration an existing disordered structure, not necessarily having the same relative concentration of component species. Since our starting system for a given simulation does not necessarily have the same density or it was generated within the required XC functional, one should make sure that there is no correlation with the initial structure. The thermal annealing cycle as presented in figure [6.16] is a reliable technique to obtain a final glass structure fully decorrelated from the starting one. Typically, to construct our amorphous model and achieve optimal statistical sampling we implemented a thermal cycle featuring 40 ps at T = 300 K, 50 ps at T = 600 K, 150 ps at T = 900 K (to allow significant diffusion) and, on cooling, 70 ps at T = 600 K and 150 ps at T = 300 K, with statistical averages taken over a final portion of the trajectory length 50 ps at T = 300 K. The lasting of each plateau could slightly change from one simulation to another, but the whole duration of the cycle ranges in between 350-450 ps. Within FPMD we have two possible ways to simulate systems under pressure, the first consists on using the constant pressure Anderson technique (see section [1.5.2.1]). In this case the initial configuration is at ambient pressure one and a target pressure is imposed to the system. To this end, the cell box is allowed to evolve isotropically to adjust the system density to the applied pressure. While the pressure is kept constant, the thermal cycle will allow the system to diffuse at each temperature. The final density of the system is computed as the average density at ambient temperature. The second method consists in generating an initial configuration at a given density (larger than the one corresponding to ambient pressure) by re-scaling the atomic coordinates and the simulation box to the desired density. The volume of the system is now kept constant during the simulation. The new configuration undergoes a geometry optimisation at T=0 K. We annealed the system at different temperatures as a part of a the thermal cycle, each temperature being controlled via a Nosé thermostat. As opposed to the first method (NPT ensemble), the pressure is computed (as 1/3 of the trace of the stress tensor) after relaxing the final structure to T= 0 K. For each density (volume) a completely independent system is generated following the same recipe, all generated systems at different volumes having only in common the same initial configuration.

V

B-Hints on the diffraction theory

For a diffraction experiment, the total neutron structure factor F(Q) is defined by:

F(Q) = n ∑ α=1 n ∑ β =1 c α c β b α b β (S αβ (Q) -1) (6.1)
where n is the number of different chemical species, α and β denote the chemical species, c α and b α are the atomic fraction and the coherent neutron scattering length of chemical species α. S αβ is the Faber-Ziman partial structure factor [1], and Q is the magnitude of the scattering vector.

Since simulation results will be compared to neutron diffraction experiment with isotopic substitution, it is worthwhile to remind how the different reciprocal and real space quantities are derived. Let us suppose that diffraction patterns are measured for three samples of glassy GeSe 2 that are identical in all aspects but different in their isotopic enrichements. If the samples are 70 Ge N Se 2 , N Ge N Se 2 and 73 Ge 76 Se 2 then the total structure factors 70 N F(Q), N N F(Q), and 73 76 F(Q) are measured, respectively. Here 70 and 73 denote glassy samples with highly enriched Ge isotopes and 76 denotes the isotopic abundance of Se. N denotes the natural isotopic abundance. The total structure factor N N F(Q) can therefore be written as:

N N F(Q) = N N A(S GeSe (Q) -1) + N N B(S GeGe (Q) -1) + N N C(S SeSe (Q) -1) (6.2)
where

N N A≡ 2c Ge c Se b N Se b N Ge , N N B≡c 2 Ge b 2 N Ge and N N C≡ c 2 Se b 2
N Se are the weighting factors (see table [6.8] for numerical values). In order to simplify the complexity associated with a single total structure factor we can form the first difference functions. For example by taking the two most contrasting structure factors In practice, from our simulations we are able to compute the Faber-Ziman partial structure factors directly in reciprocal space. Consequently we could construct all the difference function using the S αβ (Q) weighing coefficients as obtained from the experiments. This will allow to compare the same quantities from experimental and theoretical sides.

∆F Ge (Q) ≡ 70 N F(Q) - b 2 N Se b 2 76 Se 73 76 F(Q) (6.5) = F(S GeSe (Q) -1) + G(S GeGe (Q) -1) (6.6
The real space functions G(r), ∆G Ge (r) and ∆G Se (r), where r is a distance in real space, are VIII IX B-Hints on the diffraction theory obtained by using the following Fourier transform relations:

S GeGe (Q) (barns) S GeSe (Q) (barns) S SeSe (Q) (barns)
G(r) = 1 2π 2 rn 0 ∞ 0 QF(Q)M(Q) sin(Qr)dQ ∆G Se (r) = 1 2π 2 rn 0 ∞ 0 Q∆F Se (Q)M(Q) sin(Qr)dQ ∆G Ge (r) = 1 2π 2 rn 0 ∞ 0 Q∆F Ge (Q)M(Q) sin(Qr)dQ (6.9)
where M(Q) is a modification function introduced to account for the finite Q range that a diffractometer can reach (

Q max ). M(Q) is defined by M(Q) = 1 for Q ≤ Qmax , M(Q) = 0 for Q > Qmax. n 0 is the atomic number density.
Alternatively, the relevant G(r) or ∆G X (r) functions, with X=Se, Ge, Se ′ (see equation: [6.7]), could be directly obtained by replacing each S αβ (Q) by its corresponding partial pair distribution function g αβ (r) in the equations for F(Q) or ∆F X (Q), respectively.

In the diffraction theory, sets of partials structure factors exist. In what follows a reminder of the linear transformations between the different families is given: IX B-Hints on the diffraction theory X 1-From Ashcroft-Langreth to Faber-Ziman

a 11 (Q) = [S 11 (Q) -c 2 ]/c 1 (6.10) a 22 (Q) = [S 22 (Q) -c 1 ]/c 2 (6.11) a 12 (Q) = S 12 (Q)/(c 1 c 2 ) 1/2 + 1 (6.12) (6.13) 2-From Ashcroft-Langreth to Bhatia-Thornton S NN (Q) = c 1 S 11 (Q) + c 2 S 22 (Q) + 2(c 1 c 2 ) 1/2 S 12 (Q) (6.14) S CC (Q) = c 1 c 2 [c 2 S 11 (Q) + c 1 S 22 (Q) -2(c 1 c 2 ) 1/2 S 12 (Q)] (6.15) S NC (Q) = c 1 c 2 [S 11 (Q) -S 22 (Q) + (c 2 -c 1 )S 12 (Q)/(c 1 c 2 ) 1/2 ] (6.16) (6.17)
3-From Faber-Ziman to Ashcroft-Langreth

S 11 (Q) = 1 + c 1 [a 11 (Q) -1] (6.18) S 22 (Q) = 1 + c 2 [a 22 (Q) -1] (6.19) S 12 (Q) = (c 1 c 2 ) 1/2 [a 12 (Q) -1] (6.20) (6.21)
4-From Faber-Ziman to Bhatia-Thornton The PCF is given by the ratio between the density dn i j in a shell of thickness dr and the total density of the system: g ij (r) = dn ij (r) 4πr 2 drρ i (6.34) where ρ i = N i V , V is the volume of the system and N i is the number of atoms of a given species i. The position of the first peak of each g(r) i j can be taken as a measure of the bond length r i j . Furthermore, by integrating the first peak up to the first minimum, one obtains the partial coordination number n i j . XIII C-Analysis techniques and definitions XIV

S NN (Q) = c 2 1 a 11 (Q) + c 2 2 a 22 (Q) + 2c 1 c 2 a 12 (Q) (6.22) S CC (Q) = c 1 c 2 [1 + c 1 c 2 (a 11 (Q) + a 22 (Q) -2a 12 (Q))] (6.23) S NC (Q) = c 1 c 2 [c 1 (a 11 (Q) -a 12 (Q)) -c 2 (a 22 (Q) -a 12 (Q))] (6

2-Coordination number

For a given system the partial coordination numbers n i j (see Fig. [6.18]) are obtained by integrating the first peak of the corresponding g i j (r) up to the first minimum. Then, by taken into account the stoichiometry of the system one could obtain the atomic coordination numbers and the mean coordination number. The coordination number for a given species i is:

n i = n ii + n i j . (6.35)
The mean coordination number is defined by: n = c i n i + c j n j , (6.36)

where n i is the coordination number of species i and c i is its concentration.

XIV XV C-Analysis techniques and definitions

In the case of a binary system: say, GeSe 2 , we have:

n GeSe =

3-Rings

The ring statistics adopted in the present thesis work are computed via a counting algorithm based on the shortest-path criterion proposed by King and improved by Franzblau. [3,4].

Starting from an initial atom, a ring is the shortest-path between this given atom and two of its nearest neighbors.

4-Voids

We define a cavity on the basis of the generalized Voronoi-Delaunay analysis [5,6]. Atoms are considered as balls of radius r (in our case r is 50% of the average bond length for each species). This method allows the counts for multicomponent nondegenerate systems containing different atoms with different radii. Nondegenerate consider that an empty sphere cannot be tangent to more than 4 neighboring spheres.

1-First, we define Voronoi S-surfaces as the geometric locus of points equidistant from the surface of two balls. If the two balls have the same radii this surface is simply a plane, otherwise it is hyperboloid. 2-Second, we define Voronoi S-channel as the locus of points equidistant from the surface of three balls. It is also the intersection of two Voronoi S-surfaces. For three balls of equal radius the Voronoi S-chanel is a straight line. Each VS channel is characterized by a bottleneck radius R b which determines the maximum size of a test particle moving along this channel.

3-Third, we define Voronoi S-vertex as the center of an empty sphere tangent (at the same time) to four balls. In case that the four balls have the same radius we have only one vertex, balls of different radius form either more than one vertex or no vertex at all. To each vertex we associate an interstitial sphere of radius R i , tangent to the neighboring atoms. The generated network of VS vertexes and VS channels can be used to distinguish interatomic voids and resolve connected cavities. In practice, we define a void through the value of the radius R p of a probe ball (sphere) that we can insert in a given void. In addition, if a probe particle can move along an S-network bond (pass through a bottleneck (S-chanel of radius R b ), then both VS vertexes connected by this bond belong to the same cavity. The volume of voids is computed numerically, by building the VS network and taking into account also VS vertexes identified through the R i criterion as new particles with radius R i . The volume of each void is calculated numerically through a local mesh, dense enough to minimize the error below a chosen threshold. XVIII XIX C-Analysis techniques and definitions

4-Local order parameter

The local order parameter describes the atomistic geometrical arrangement within a network. Representative values are q = 1 for the ideal tetrahedral geometry and q = 0 for the six-coordinated octahedral sites. Values in between q=0 and q=1 refer to a given defective octahedral configurations. The q-parameter is defined by: q ≡ 1 -3 8 ∑ k>i [ 1 3 + cosθ i jk ] 2 where θ i jk is the angle formed between a central atom j and its neighboring atoms i and k [7,8]. Corner-sharing configuration: This is when two Ge atoms share one Se atom (Ge 3 , Ge 4 , Ge 5 , and Ge 6 in fig. [6.22]).

Homopolar bonds: This is when two atoms of the same chemical species are directly connected in a nearest-neighbor configuration to form a bond (Ge 6 -Ge 7 and Se 1 -Se 2 in fig. [6.22]). When more than two atoms of the same chemical species connect to each other without exceeding a coordination 2 (except the for the first and the last atoms) they form a homopolar chain (Se 4 -Se 5 -Se 6 ).

n-fold X atom: When a given atom (in the above figure) makes n bonds (Ge 8 : 2-fold Ge, Ge 10 : 3-fold Ge, Ge 9 : 4-fold Ge, Ge 11 : 5-fold Ge, Se 5 : 3-fold Ge, Ge 3 : 4-fold Ge).

XX

E-Parrinello-Rahman barostat

The Parrinello-Rahman method [9][10][11][12] is the natural extension of the Anderson barostat to the non isotropic external stress case. The substantial intake of Parrinello and Rahman is the simultaneously control of the shape and the volume of the simulation box that change as a consequence of the external pressure. The change of the lattice vectors corresponds to introducing them as a new dynamical variable to the Lagrangian. In practice, the real space position of a given particle R I is expressed in term of primitive Bravais lattice vectors a, b and c of the simulation box.

For simplicity it is better to work with scaled coordinates that do not depend on the the real cell volume. Thus we define the scaled variable s as: r = sh (6.50) where h is the 3 × 3 matrix representation of the Bravais lattice vectors h = (a, b, c). In this framework the original volume of the cell is given by Ω = det(h) and the conversion of distances from real to scaled coordinates follows:

(r i -r j ) 2 = (s i -s j ) t G(s i -s j ) (6.51)

Where G is the conversion metric tensor defined as: S I is the scaled coordinate that defines the I th particle position on the unit cube Ω unit = 1. The presence of the electronic wave function in the CP Lagrangian makes the introduction of the variation on the cell shape delicate to treat. One way to get around this problem is to express the orbitals in term of the scaled variable s. In this case the original wave functions φ i is such that:

G = h T h =   
φ i (r) = 1 √ Ω φ i (h -1 r) = 1 √ Ω φ i (s) (6.54)
The new wave function φ i (s) is independent on the Bravais lattice matrix h and the normaliza- tion relation is invariant under scale transformation: As a consequence, the electronic charge density could also be rewritten in term of the scaled scaled variable s:

n(r) = 1 Ω n(s) (6.56)
This equation demonstrates that the scaled wave functions and charge density are independent on the change on the dynamical variable associated to the cell degree of freedom. Also from Eq.[6.54] and Eq. [6.56] one could read that the only direct effect of the change on the cell shape on both the wave function and the electronic density is to change them so as to be compatible with the boundary conditions. In this manner, all variables could be written in a scaled cell-independent fashion, and one could write the extended version on the CP Lagrangian (Eq. [1.43]):

L = 1 2
µ ψ * i (s) ψi (s)ds -E KS ({ψ i }, {hS I }) + ∑ i j λ i j ( µψ * i (s)ψ j (s)ds -δ i j )

+ ∑ I 1 2 M I ṠT I G ṠI + 1 2
W Tr ḣT ḣ -pΩ (6.57)

Notice here that the all the integrals in Eq. [6.57] are taken over the scaled cell with unitary volume Ω unit , p defines the external hydrostatic applied pressure, W is the fictitious mass parameters that we associate to the time scale evolution of the cell motion h. In case of rigid cell, ḣ → 0, the PR Lagrangian reduces to the initial CP one (Eq. The total internal stress tensor reads:

Π tot us = 1 Ω ∑ I M I ṠT I G ṠI + Π us (6.61)
which is the sum of the thermal contribution of the nuclear motion at a temperature T and the electronic internal stress tensor Π defined as:

Π uv = - 1 Ω ∑ s δ E tot δ h us h T sv (6.62)
The average internal pressure, calculated as (1/3)Tr(Π tot ) , should be equal to the external pressure p, since from Eq.[6.60] one notice that the total internal stress Π tot should be dynamically proportional to Pδ . This is ensured due to the friction coefficient ∝ Ġ in Eq. [6.58], recalling the frictional feedback mechanism showed in the case of Nosé-Hoover thermostat.

XXIII

F-GPW basis set

Density functional theory (DFT) as introduced in Part I of the manuscript is an efficient method to perform electronic structure calculations. The predictive-power of the DFT give access to a large range of properties for a wealth of systems ranging from solid state physics to biochemistry. Nevertheless, the size of the system can be a stumbling block still preventing applications in specific cases requiring a minimal number of atoms to allow a meaningful description of certain properties. In its Kohn-Sham expression, most DFT implementations are highly non-linear as a function of the system size, leading to an expansive computational cost. Such high computational cost is dominated by the orthogonalization of the wave functions and the calculation of the Hartree energy [13]. An alternative technique to circumvent these problems has been proposed by Lippert, Hutter and Parrinello [14]. It consists on the use of a mixed hybrid Gaussian and plane wave basis set (GPW method, originally first implemented in the CP2k code). Indeed, one of the ternary systems studied in this theses was generated using the GPW method within the CP2k suite of program. Since DFT calculations as implemented in CP2k are (mostly) based on GPW basis set, it is worthwhile to give some details on this method. GPW is a hybrid scheme using an atom-centered Gaussian-type basis set to expand and describe the wave functions and an auxiliary plane wave basis set to describe the density. Contrary to the plane waves, the use of a localised Gaussian reduces considerably the size of the basis set and lead to a considerable reduction on the computational cost. The plane wave descriptions to the density allows to compute the Hartree energy by means of Fast Fourier transforms. As a consequence of this duality, the computational cost needed to obtain the Kohn-Sham matrix scales linearly with the system size.

In short, we detail hereafter the expression for the density, that is based on the sum of two distinct representations. We give the expressin for the total density as well. The electron density n(r) is expanded on an atom centered, contracted Gaussian functions: where Ω is the unit cell volume, G are the reciprocal space vectors and ñ(G) are the expansion coefficients ñ(G). ∑ A n A (r) represents the hard contributions from the atomic cores expended in Gaussian basis set:

n A (r) = ∑ µν∈A P µν A g µ (r)g * ν (r) (6.66)
The ñA (r) represents a smooth local term used to compensate for the overlap between the soft and hard densities in the atomic core regions, it is also expanded in Gaussian basis set: ñA (r) = ∑ µν∈A Pµν A g µ (r)g * ν (r) (6.67)

Figure [6.23] shows an illustration of the wave function expansion with the GPW method. Near the the atom A the density is given by the density of the atom A expanded in the localised Gaussian basis set (n(r) = n A (r)). However, far from the atom the density is equal to the one expanded on a plane wave basis set. The way the density is expanded in the GPW method allow the efficient and rapid conversion between n(r), ñ(r) and ñ(G) by mean of FFT and a reliable mapping procedure. Within this representation the KS DFT energy (Eq.1. This part and the related extensions (pseudo potentials, XC functionals, ..) is the DFT part of the CP2k code [15] implemented in the QUICKSTEP scheme to efficiently solve the Poisson equation. A detailed description of this scheme could be found in ref. [16].

XXVII

G-HSE06 calculations

In the attempt to improve the level of agreement between the FPMD and the ND at high temperature, we performed additional simulation using the generalization of the PBE0 [17] functional proposed by the group of Scuseria [18] (called HSE06), in which the PBE exchange functional [19] is complemented by the Hartree exact exchange with a properly calibrated mixing coefficient. The aim is to ascertain whether an increased accuracy in the description of the bonding character is instrumental in better reproducing the experimental outcome at high pressure. This recipe was applied for a single system at V /V 0 =0.654 (ρ/ρ 0 = 1.529). We have performed a thermal cycle featuring 1.2ps @300K then 3.42ps@1100 then 5.80ps@900K and finally 6.58ps@300K. The results were obtained from analysing the relaxed trajectory at 300K during the cooling cycle. The obtained pair correlation function is compared to the one obtained from the BLYP simulation at the same reduced density (the same system). The following figure illustrates the comparison: XXIX G-HSE06 calculations

XXX

The shape of the partial pair correlation function computed with HSE06 and BLYP show a some differences. In particular the short range is affected. The first peak in g HSE06 GeGe (r) is shifted toward higher distances and no well defined minimum could be distinguished. Looking at g HSE06 SeSe (r), the first peak is higher than the BLYP one resulting on a higher fraction homopolar of bonds and a larger partial coordination number. In contrast, the first peak intensity in g HSE06 GeSe (r) decreases when compared to the BLYP result. In addition the main peak position shifts leftward. HSE06 mean coordination gives number comparable to that found from BLYP model while the HSE06 mean bond length is by far larger than its BLYP homologue. This large bond length is related to the treatment of exchange part in the hybrid XC functional. In conclusion, the account of a higher level XC recipe does not improve by itself the level of agreement between the FPMD and the neutron diffraction experiment.

XXX

H-Extra Phd contributions

The work performed during this Phd project led to the publication of papers in international recognized scientific journals. Beside this work, I contributed to other research projects within our team. The papers related to the Phd and the extra-Phd contributions (published and in preparation) are sorted below. For the extra-Phd contributions, abstracts are recalled. Further details could be found on the online version of the papers. The structural properties of glassy Ge 2 Se 3 were studied in the framework of first-principles molecular dynamics by using the Becke-Lee-Yang-Parr scheme for the treatment of the exchange-correlation functional in density functional theory. Our results for the total neutron structure factor and the total pair distribution function are in very good agreement with the experimental results. When compared to the structural description obtained for liquid Ge 2 Se 3 , glassy Ge 2 Se 3 is found to be characterized by a larger percentage of fourfold coordinated Ge atoms and a lower number of miscoordinations. However, Ge -Ge homopolar bonds inevitably occur due to the lack of Se atoms available, at this concentration, to form GeSe 4 tetrahedra. Focusing on the family of glasses Ge x Se 1-x , the present results allow a comparison to be carried out in reciprocal and real space among three prototypical glassy structures. The first was obtained at the stoichiometric composition (glassy GeSe 2 ), the second at a Se-rich composition (glassy GeSe 4 ) and the third at a Ge-rich composition (glassy Ge 2 Se 3 ). All networks are consistent with the 8 -N rule, in particular, glassy GeSe 4 , which exhibits the highest degree of chemical order. The electronic structure of glassy Ge 2 Se 3 has been characterized by using the Wannier localized orbital formalism. The analysis of the Ge environment shows the presence of dangling, ionocovalent Ge-Se, and covalent bonds, the latter related to Ge-Ge connections.

❅ The structure of liquid GeSe revisited: A first principles molecular dynamics study. J. Chem. Phys. 138, 174505 (2013); Sébastien Le Roux, Assil Bouzid, Mauro Boero, and Carlo Massobrio.

Early first-principles molecular dynamics results on liquid GeSe were characterized by shortcomings in the description of Ge -Ge (and to a lesser extent Se -Se) short range correlations. In that case the exchange-correlation functional adopted was the one devised by Perdew and Wang (PW91). In the search of improvements in the atomic-scale modelling of this liquid, we have produced new sets of data by employing two different schemes for the exchange-correlation part within the density functional theory approach. The two functionals selected are those proposed by Becke, Lee, Yang, and Parr (BLYP) and by Perdew, Burke, and Ernzerhof (PBE). The PBE results turned out to be quite similar to the PW91 ones. The XXXII XXXIII H-Extra Phd contributions BLYP results feature instead a better account of the Ge-Ge first shell of neighbors, correctly exhibiting two clear maxima separated by a deep minimum. Due to the increase in the number of the tetrahedral structural units, the atomic mobility of Ge and Se atoms in the network is reduced with respect to the PW91 case. This brings the diffusion coefficients of the two species down to values close to those of liquid Ge 2 Se 3 and liquid GeSe 2 . The structure of glassy GeS 2 is studied in the framework of density functional theory, by using a fully self-consistent first-principles molecular dynamics (FPMD) scheme. A comparative analysis is performed with previous molecular dynamics data obtained within the Harris functional (HFMD) total energy approach. The calculated total neutron structure factor exhibits an unprecedented agreement with the experimental counterpart. In particular, the height of the first sharp diffraction peak (FSDP) improves considerably upon the HFMD results. Both the Ge and the S subnetworks are affected by a consistent number of miscoordinations, coexisting with the main tetrahedral structural motif. Glassy GeS 2 features a short-range order quite similar to the one found in glassy GeSe 2 , a notable exception being the larger number of edge-sharing connections. An electronic structure localization analysis, based on the Wannier functions formalism, provides evidence of a more enhanced ionic character in glassy GeS 2 when compared to glassy GeSe 2 . First-principles calculations within the framework of the density functional theory are used to build realistic models for the surface of glassy GeS 2 (a -GeS 2 ) at 0 K and 300 K. This allows comparison between the structural and electronic properties of surface and bulk a -GeS 2 . Thanks to surface bond rearrangements (thermally-induced), both Ge and S atoms increase their average coordination at the free surface. Although the a -GeS 2 surface recovers the main tetrahedral structural motif of bulk a -GeS 2 , the a -GeS 2 surface shows a slightly lower chemical order. This shift with respect to the ideal stoichiometry is balanced by a greater content of overcoordinated S atoms and a greater content of three-fold coordinated Ge atoms with respect XXXIII H-Extra Phd contributions XXXIV to the percentage of two-fold coordinated Ge. This effect is fostered when the a -GeS 2 surface is obtained at 0 K. Maximally-localized Wannier functions (WF) are employed to characterize the nature of chemical bonding in the structural units of the a -GeS 2 surface. We compare the ability of several charge derivation methods to capture the atomic charge variations induced by a coordination change. Our estimate for the charges allows exploiting the first-principles results as a database to construct a reliable interatomic forcefield. XXXIV structure des chalcogénures binaires et ternaires sous différentes conditions thermodynamiques. Des techniques de modélisation numérique ab-initio ont été employées. En particulier, nous avons utilisé la dynamique moléculaire par premiers p i ipes selo l'app o he de Ca et Pa i ello ai si ua sa ve sio dite "deu i e g atio ". La p e i e pa tie est o sa e à l' tude des chalcogénures binaires, notamment les verres GeSe 2 et GeSe 4 sous p essio ai si u'à l' tude des effets des fo es de va de Waals et des fo tio elles d' ha ge et o latio DFT su la st u tu e de l'a o phe GeTe 4 . Da s la deu i e pa tie, l'i t t a t po t à l' tude d'u atériau à changement de phase récemment proposé par les expérimentateurs comme un bon candidat pour le stockage de données, le ternaire Ga 4 Sb 6 Te 3 .

❅ First-principles molecular dynamics study of glassy

❅ Surface of glassy

Abstract :

This thesis reflects efforts toward an accurate understanding of the atomic scale structure of chalcogenide glasses. These compounds have an impact on electronics, optoelectronics and memory devices. I resorted to the most advanced first-principles molecular dynamics simulations such as the standard Car-Parrinello method as well as its second generation version. In the first part of this thesis we provide a detailed study of the topological changes undergone under pressure by glassy GeSe 2 and by glassy GeSe 4 . Structural transition and bonding features are described and compared to the results of neutron and X-ray diffraction experiments. Furthermore, in the case of glassy GeTe 4 we demonstrated that the inclusion of van der Waals forces leads to substantial improvements in the description of the structure. In the second part of this thesis, we established the atomic-scale organization of a promising candidate for phase change memory applications, glassy Ge 4 Sb 6 Te 3 .

  the thermal frequencies of the ionic and the electronic sub-system, as well as theNosé thermostats for both sub-systems in CPMD. . . . . . . . . . . . . . . . . . . . . 1.10 (a): Ambient pressure,time evolution of the fictitious Car-Parrinello electron kinetic energy (black line) and the ions kinetic energy (red line) in a typical CPMD simulation at different temperature. (b): High pressure (9.8GPa), time evolution of the fictitious Car-Parrinello electron kinetic energy (black line) and the ions kinetic energy (red line) in a typical CPMD simulation at different temperature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Total energy conservation for different barostat fictitious mass. . . . . . . . . . . . . . . . . . . . 1.12 Time evolution of the cell side length for different barostat fictitious mass. . . . . . . . . . . . . . 1.13 Time evolution of the diagonal and the off-diagonal compounds of the stress tensor during the MD simulation at T= 300K. The annealing region corresponds to a T=0 K simulation. . . . . . . . . .

3 . 4 70 GeSe 2 .

 34702 and Co. ((black) • with vertical error bars) are compared to the results obtained from first-principles molecular dynamics in the NVT present work ((red) ), NPT present work results ((magenta) ×) and in the work by Durandurdu and Drabold ((blue) ♦). The computed and the measured data are fitted to a second-order Birch-Murnaghan equation of state (solid blue and solid black curves respectively). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii The pressure dependence of the total structure factor (a) F(k) for N N GeSe 2 and (b) 70 N F(k) for N The dark solid (black) curves show spline fit of the the data sets measured using (a) the PEARL diffractometer, except for the range k < 1.55 Å where they represent fitted Lorentzian functions because this region was not accessible in the diffraction experiments. (b) using the D4C diffractometer. The solid light (green) curves show the molecular dynamics results for densities comparable to those of the measured samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 The pressure dependence of (a) ∆F Ge (k) and (b) ∆F Se (k). The data sets, measured using D4 (points with vertical error bars), are compared to the FPMD results at similar densities (light (green) curves). The dark (black) curves are the back Fourier transforms of the measured r-space difference functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.

  4 , where Ge atoms are dark (blue) and Se atoms are light (green). Se atoms along a connection path between two Ge atoms are labelled as AA, those between one Ge atom and one Se atom are labelled as AB, and those between two Se atoms are labelled as BB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.23 (a): The pressure dependence of the fractions of 2-fold Se broken down into the contributions of Se of type AA, AB and BB. (b) and (c): The pressure dependence of the fractions of l-fold Ge atoms (l=3, 4 or 5) and l-fold Se atoms (l= 1, 2 or 3). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.24 Bond angle distribution as a function of pressure for Ge -Se -Ge, Se -Ge -Se, Ge -Se -Se and Se -Se -Se. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 3.25 Atomistic configurations of glassy GeSe 4 system at ambient and high pressure. The blue arrow indicates the Se 3 connection, the red one the Ge 5 connection and the green one indicates the Ge -Ge wrong bond. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Top figure: amorphous GeTe 4 total structure factor. Bottom figure amorphous GeTe 4 total pair correlation function.

  ) and model C (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 The partial pair correlation functions of amorphous Ga 4 Sb 6 Te 3 . Model A (green line), model B (blue line) and model C (red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Bond angle distributions: Sb-Ga-Sb, Te-Ga-Te, Ga-Sb-Ga, Te-Sb-Te, Ga-Te-Ga and Sb-Te-Sb. . . . 6.8 Local order parameter distribution for the 3-and 4-fold Ga, Sb or Te atoms. Model A (black lines), model B (red lines) and model C (blue lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Rings distributions for the amorphous phase. Model A (red), model B (green) and model C (blue) . 6.10 (upper plot) Snapshots of C model. Only Sb atoms (yellow) and Te atoms (green) are shown. (bottum plot) Snapshot of the 97-atom Sb-cluster in the 299-atoms model of Ga 4 Sb 6 Te 3 . The atoms of the cluster are shown as black spheres while all the other atoms are represented by gray (red), black and light gray (cyan) lines for Ga, Sb and Te respectively. The shape of the cluster is depicted by the (blue) transparent surfaces of the Voronoi polyhedra built around the Sb atoms of the cluster. Four 299-atoms supercells are shown but the periodic boundary condi-tions are not applied to the atoms of the Sb cluster to highlight its shape and connectivity . . . . . . . . . . . . . . . . . xi 6.11 Distribution of the volume of cavities (left plot) and distribution of Sb clusters (right plot) in model C. The cluster size is indicated as the number of atoms per cluster. . . . . . . . . . . . . . . . . . 178 6.

  Two-dimensional illustration of the Voronoi regions in a system of balls of different sizes. The edges of the Voronoi polyhedra are shown by thick lines. The dashed lines depict the edges of the conventional Voronoi polyhedra of the system as computed for the centers of the balls. Right figure: Different types of configurations of four balls in 3D (left) and three disks in 2D (right). (a) Singlet set: one inscribed sphere; (b) Double set: two inscribed spheres; (c) Nonsimpliciable: no inscribed spheres. This figure contains two figures shown as separate in ref. [5]. . . . . . . . . XVIII 6.21 Snapshots of tetrahedral, octahedral and defective octahedral configurations. The value of the local order parameter is provided. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIX 6.22 A snapshots of a three-dimensional view of structural subunits of amorphous GeSe 2 . . . . . . . . XX 6.23 Wave function expansion within the GPW approximation. . . . . . . . . . . . . . . . . . . . . . . XXVII 6.24 Partial pair correlation functions as computed with HSE06 and BLYP functionals. . . . . . . . . XXIX xii List of Tables 1.1 Energy conservation in a.u./ps for Car-Parrinello (CPMD) and Born-Oppenheimer (BOMD)

  = 10.0(1), b 73 Ge = 5.09(9) and b 76 Se = 12.2(1) fm [2]. The atomic fraction of Ge and Se are c Ge = 1/3 and c Se = 2/3 respectively. . . . . . . . . . . . . . . IX xiv First of all I would like to thank my supervisor Carlo Massobrio, I have learned from him an enormous amount not only about research and science, but also about communication, life, politics and the need of balance in all things. You offered me full support; the door of your office was always open for me. Thank you for your patience, for encouraging and pushing me up. Thank you for the trust you gave me and for striking just the right balance between advising me and enabling me to explore interesting side avenues.

Figure 2

 2 Figure 2Cycle thermique: l'évolution de la température au cours de la simulation.

Figure 4

 4 Figure 4 Configurations atomistiques du verre GeSe 2 prises à différentes pressions. Couleur violette: Atomes de Ge et couleur jaune: atomes de Se.

Figure 5

 5 Figure 5 Fonctions de corrélation de paires partielles pour l'amorphe GeSe 4 sous pressions.Du haut jusqu'en bas: g SeSe (r), g GeSe (r), g GeGe (r) et fonction de corrélation radiale totale g tot (r).

Figure 6

 6 Figure 6 Configurations atomistiques du verre GeSe 4 à pression ambiante et à 13 GPa. La flèche bleue indique une configuration de type Se 3 , la rouge une configuration de type Ge 5 et la verte indique un défaut de type Ge -Ge.

Figure 7

 7 Figure 7 Fonctions de corrélation de paires partielles pour l'amorphe GeTe 4 . Du haut jusqu'en bas: g GeGe (r), g GeTe (r) et g TeTe (r) De gauche à droite: modèle généré avec PBE et modèle généré avec BLYP. Simulations avec les corrections de vdW (lignes rouges) et les simulations sans les corrections de vdW (lignes noires).

Figure 8

 8 Figure 8 La distribution de paramètre d'ordre local pour le modéle BLYP+vdW (couleur rouge) et pour le modéle PBE+vdw (couleur noire). Les flèches indiquent les positions qui correspondent à des configurations octaédriques défectueuse. Des photos du modéle PBE+vdw sont également fournis.

Figure 9 L

 9 Figure 9 L'évolution de la température au cours de la trempe du modèle A (couleur verte en pointillés), du modèle B (couleur bleue) et du modèle C (couleur rouge).

Figure 10

 10 Figure 10 Configurations atomistiques du verre Ga 4 Sb 6 Te 3 avec les principaux motifs structuraux.

13 1. 1 .

 131 Modeling materials by classical molecular dynamics : advantages and pitfalls Euler-Lagrange equation: d dt

Figure 1 . 1

 11 Figure 1.1 General scheme of a computer code for MD

  E xc have been proposed to include functional with different complexity levels. Perdew proposed a way to classify the E xc functional [22], Fig.[Fig 1.2] illustrates the classification.

Figure 1 . 2

 12 Figure 1.2 Schematic diagram of "Jacob's ladder" of exchange correlation functionals proposed by J. P. Perdew.

Figure 1 . 3

 13 Figure 1.3 Schematic illustration of all-electron (dashed lines) and pseudopotential (solid lines) and their corresponding wave functions. The radius at which the allelectron and pseudo-electron values match is designated as r c .

31 1. 4 .

 314 Car-Parrinello molecular dynamics parameter µ increases. Similarly the highest electron frequency is given by:

Figure 1 . 4

 14 Figure 1.4 Heat transfer as a function of the prefixed average fictitious kinetic energy E kinof the electronic wave functions for solid aluminum at the melting point. The arrow indicates the kinetic energy required for adiabatic motion of the electrons according to eq.[1.72]. Taken from ref.[75]

Figure 1 . 6

 16 Figure 1.6 Car-Parrinello conserved total energy variation for different temperatures. Tests performed on a GeSe 2 system.

Figure 1 . 7

 17 Figure 1.7 The simulation time per molecular dynamic step variation as a function of number of CPUs used. (a): a system of Ga 4 Sb 6 Te 3 with 117 atoms has been simulated on IBM x3750 M4 (Ada-IDRIS). (b) results of Ga 4 Sb 6 Te 3 simulation with 299 atoms on SGI ICE machine (Jade-CINES). (c) results of Ga 4 Sb 6 Te 3 simulation with 299 atom on IBM Blue Gene / Q (Turing-IDRIS)

Figure 1 . 8

 18 Figure 1.8 Temperature variation and electronic fictitious kinetic energy as imposed by the velocity rescaling method.

Figure 1 . 9

 19 Figure 1.9 Schematic representation of the thermal frequencies of the ionic and the electronic sub-system, as well as the Nosé thermostats for both sub-systems in CPMD.

Figure 1 .

 1 Figure 1.10 (a): Ambient pressure,time evolution of the fictitious Car-Parrinello electron kinetic energy (black line) and the ions kinetic energy (red line) in a typical CPMD simulation at different temperature. (b): High pressure (9.8GPa), time evolution of the fictitious Car-Parrinello electron kinetic energy (black line) and the ions kinetic energy (red line) in a typical CPMD simulation at different temperature.

Figure 1 . 11

 111 Figure 1.11 Total energy conservation for different barostat fictitious mass.

  good stability is recovered. In the following figure we display the cell parameter variation for different parameters W . For W = 10 6 a.u. the cell dynamics oscillates very rapidly around

Figure 1 . 12

 112 Figure 1.12 Time evolution of the cell side length for different barostat fictitious mass.

Figure 1 . 13

 113 Figure 1.13 Time evolution of the diagonal and the off-diagonal compounds of the stress tensor during the MD simulation at T= 300K. The annealing region corresponds to a T=0 K simulation.

Figure 2 . 1

 21 Figure 2.1 Schematic of the three elastic phases observed in network glasses as a function of increasing connectivity. In the cases of truly random networks the intermediate phase collapses yielding an elastic phase transition from a floppy to a stressed rigid. Taken from [1].

  [2.2]) and a non negligible fraction of wrong bonds has been reported.

Figure 2 . 2

 22 Figure 2.2 Upper curve, static structure factor for 63-atom model. The lower curve and points are experimental results from ref [6]. Figure taken from [2].

2 Figure 2 . 3

 223 Figure 2.3 Faber-Ziman partial structure factors for liquid GeSe 2 : GGA calculations (solid line), LDA calculations (dotted line), and experiment (dots with error bars) Ref.[12] . S LDA GeGe (k), S LDA GeSe (k), and S LDA SeSe (k) have been shifted down by 2, 1, and 1, respectively. From ref[11]

Figure 2 . 4

 24 Figure 2.4 Contour plots of valence electronic charge density for a selected Se-Ge-Se trimer within liquid GeSe 2 in the LDA (upper left) and GGA (upper right). A difference plot (GGA -LDA) is also shown (lower). The scales are in atomic units. FromRef.[13]

GeSe 2 Figure 2 . 5

 225 Figure 2.5 Calculated neutron structure factor for liquid GeSe 2 , obtained within GGA (solid line) and LDA (dots), compared to experimental data (circles) ref.[14]. For clarity, LDA curve is displaced downward by 0.4. From ref.[15].

  [2.6].

Figure 2 . 6

 26 Figure 2.6 Calculated neutron structure factor for liquid GeSe 2 (black lines) compared to experimental data (circles). From ref. [9].

67 2. 2 .

 672 GeSe 2 Overall, this work provided a remarkable improvement of the microscopic description of the structure of liquid GeSe 2 and a good agreement with the diffraction data. However, it did not contribute to explain the physical origin of the FSDP in the S cc (k), which remains (at the level and at the moment of this specific publication) still an open issue.

  [2.7]).

Figure 2 . 7

 27 Figure 2.7 Total neutron structure factors for liquid GeSe2 at T = 1050 K and T = 1373 K.Open circles: experimental results, from[12] (T = 1050 K) and[20] (T = 1373 K). Solid line: theory,[9] (T = 1050 K) and[19] (T = 1373 K). From Ref.[19]

Figure 2 . 9

 29 Figure 2.9 Partial pair-correlation functions for liquid GeSe 2 : black line: present BLYP results, red line: PW results of Ref.[11] , open circles: experimental results of Ref.[12]. From Ref.[26]

Figure 2 . 10

 210 Figure 2.10 Representative subset of Ge and Se atoms in amorphous GeSe 4 where Ge atoms are dark (blue) and Se atoms are light (green). Se atoms along a connection path between two Ge atoms are labelled as AA, those between one Ge atom and one Se atom are labelled as AB, and those between two Se atoms are labelled as BB.

Figure 2 . 11

 211 Figure 2.11 The calculated partial pair correlation functions for glassy GeSe 4 . The label Se-Se indicates g SeSe (r), the label Ge-Se indicates g GeSe (r) and the label Ge-Ge indicates g GeGe (r). The circle symbols curve corresponds to PW, the other one to PBE. A set of representative error bars are given in the g GeGe (r) case. Figure taken from Ref.[32]

Figure 2 .Figure 2 . 13

 2213 Figure 2.12 Pressure-volume equation of state for amorphous GeSe 2 . Experimental compression in a hydrostatic medium, a 4:1 methanol:ethanol mixture (solid circles); experimental decompression from 8.55 GPa (open triangles); simulation data open squares Ref.[41] ; third-order isothermal Birch-Murnaghan equation of state fit to the experimental compression data (solid curve). From Ref.[40]

  [2.14]).

Figure 2 . 14

 214 Figure 2.14 Variation of the average coordination number of the 1st peak in g(r) glasses as a function of pressure. The open circles are glassy GeSe 4 data obtained from integrating the plotted function between the two minima either side of the first peak.The filled circles are GeSe 2 data from Ref.[40]. The blue triangles are ambient data from Ref.[44]. Figure from Ref.[43]

Figure 3 . 1

 31 Figure 3.1Partial pair correlation functions for glassy GeSe 2 . The experimental results of Ref.[5] (solid curve) are compared with the results of Ref.[1] with ρ a (dashed red curve) and with the new simulation at the experimental density ρ exp (dotted blue curve with stars) and the present relaxed structure results at ρ b (NVT ensemble, dashed-dotted green curve).

Figure 3 . 2

 32 Figure 3.2 The Bhatia-Thornton concentration-concentration partial structure factor S cc (k)for glassy GeSe 2 . Experimental results of Ref.[5] (black solid line), simulation results with ρ a Ref.[1](red dashed curve), simulation results with ρ exp (blue dotted curve and simulation results of the calibrated system at ρ b (green dashed curve).

Figure [ 3

 3 .4-a] and[3.4-b] display direct comparison between the computed and measured total structure factors for the70 Ge N Se compound and the natural compound respectively.

Figure 3 . 4

 34 Figure 3.4 The pressure dependence of the total structure factor (a) F(k) for N N GeSe 2 and (b) 70 N F(k) for N 70 GeSe 2 . The dark solid (black) curves show spline fit of the the data sets measured using (a) the PEARL diffractometer, except for the range k < 1.55 Å where they represent fitted Lorentzian functions because this region was not accessible in the diffraction experiments. (b) using the D4C diffractometer. The solid light (green) curves show the molecular dynamics results for densities comparable to those of the measured samples.

Figure 3 . 5 1 .

 351 Figure 3.5 The pressure dependence of (a) ∆F Ge (k) and (b) ∆F Se (k). The data sets, measured using D4 (points with vertical error bars), are compared to the FPMD results at similar densities (light (green) curves). The dark (black) curves are the back Fourier transforms of the measured r-space difference functions.

Figure 3 .

 3 Figure 3.6 The G(r) functions from experiment (solid dark (black) curves) were obtained byFourier transforming the measured F (k) functions shown in Fig.[3.4] via the application of a Lorch modification function. The broken dark (black) curves show the extent of the unphysical oscillations at r-values smaller than the distance of closest approach between the centers of two atoms. The solid light (green) curvesshow the Fourier transforms of the simulated functions shown in Fig.[3.4] after applying the same maximum cutoff value k max = 19.55 Å -1 as used for the neutron diffraction data.

Figure 3 . 7

 37 Figure 3.7 The G(r) functions from experiment (solid dark (black) curves) were obtained by Fourier transforming the measured F (k) functions shown in Fig.[3.4] The solid light (green) and dashed line (blue) curves show the Fourier transforms of the simulated functions of simulated systems with and without thermal cycle respectively.

Figure 3 . 8

 38 Figure 3.8 The pressure dependence of the partial pair distribution functions from FPMD simulations. From the top to the bottom: the partial radial distribution function g SeSe (r), g GeSe (r) and g GeGe (r).

  101

Figure 3 . 10

 310 Figure 3.10 The pressure dependence of the fractions of Ge CS atoms and ES Ge atoms isgiven, together with a breakdown of the latter into its contributions from Ge 1 , Ge 2 and Ge 3&4 units. Left plot: NVT simulation results, right plot: NPT simulation results.

  Figure 3.11The pressure dependence of the intra-polyhedral Se-Ge-Se and the interpolyhedral Ge-Se-Ge bond angle distributions. NVT simulations (left) and NPT simulations (right).

Figure 3 . 12

 312 Figure 3.12The ρ/ρ 0 dependence of the intra-polyhedral Se-Ge-Se ((red) with vertical er- ror bars) and inter-polyhedral Ge-Se-Ge ((black) with vertical error bars) bond angles as estimated from the measured Ge-Se, Ge-Ge and Se-Se distances. The experimental results are compared to the mean values θ SeGeSe (broken (red) curve with symbols) and θ GeSeGe (broken (black) curve with symbols) taken from the FPMD simulations.

107 3 . 1 .

 31 Amorphous germanium diselenide under pressure of each n-fold Ge atom show that the peaks in the range 0.2 to 0.6 are mainly due to 5-fold Ge either in trigonal bipyramid or square based pyramids. At ρ/ρ 0 = 1.639 the network features a distorted octahedral configuration coexisting with defective GeSe 5 trigonal bipyramids and square based pyramids.

Figure 3 . 14

 314 Figure 3.14 Pressure dependence of the n-membered rings

Figure 3 . 15

 315 Figure 3.15The pressure dependence is given for the fractions of n-fold coordinated Ge and Se atoms, along with the fractions of these n-fold coordinated Ge and Se atoms that contain homopolar bonding defects. The error bars (usually smaller than the symbol size) were calculated according to Ref.[16]. Left plot: NVT simulation results, right plot: NPT simulation results.

109 3 . 1 .

 31 Amorphous germanium diselenide under pressure coordination number and the occurrence of the pseudo cubic network.

24 Figure 3 . 16

 24316 Figure 3.16 Atomistic configurations taken from GeSe 2 glass at different pressures. Ge atoms are dark (purple) and Se atoms are light (yellow). Bonds are drawn when two atoms are separated by a cutoff distance given by the position of the first minimum in g GeSe (r).

Figure 3 . 19

 319 Figure 3.19The pressure evolution of the partial Faber-Ziman structure factor from from Fourier transform of the partial real space counterpart presented on figure.[3.20]. From the top to the bottom S FZ SeSe (k), S FZ GeSe (k) and S FZ SeSe (k).

Figure 3 . 20

 320 Figure 3.20Pair correlation function of amorphous GeSe 4 under pressures. from top to bottom g SeSe (r), g GeSe (r), g GeGe (r) and total radial correlation function g tot (r).

Figure 3 .

 3 Figure 3.22 A representative subset of Ge and Se atoms, in amorphous GeSe 4 , where Ge atoms are dark (blue) and Se atoms are light (green). Se atoms along a connection path between two Ge atoms are labelled as AA, those between one Ge atom and one Se atom are labelled as AB, and those between two Se atoms are labelled as BB.

Chapter 3 .

 3 Binary

Figure 3 .

 3 Figure 3.23 (a): The pressure dependence of the fractions of 2-fold Se broken down into the contributions of Se of type AA, AB and BB. (b) and (c): The pressure dependence of the fractions of l-fold Ge atoms (l=3, 4 or 5) and l-fold Se atoms (l= 1, 2 or 3).

GeSe 2 .Figure 3 . 24

 2324 Figure 3.24 Bond angle distribution as a function of pressure for Ge -Se -Ge, Se -Ge -Se, Ge -Se -Se and Se -Se -Se.

Figure 3 . 25

 325 Figure 3.25 Atomistic configurations of glassy GeSe 4 system at ambient and high pressure. The blue arrow indicates the Se 3 connection, the red one the Ge 5 connection and the green one indicates the Ge -Ge wrong bond.

  127

Figure 4 . 1

 41 Figure 4.1 Top figure: amorphous GeTe 4 total structure factor. Bottom figure amorphousGeTe 4 total pair correlation function. Experimental results from ref.[5] (black dots), BLYP results (blue solid lines), BLYP+vdW results (red solid lines), PBE results (purple dashed lines) and PBE+vdW results (green solid lines). PBE and PBE+vdW was shifted to y+2 (for the s(k)) and to y+2.5 (for the g(r)).

Fig.[ 4 . 2 ] 4 . 3 .Figure 4 . 2

 424342 Fig.[4.2] displays the calculated partial correlation functions for the amorphous GeTe 4 treated

Chapter 4 .

 4 Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4

Figure 4 . 3

 43 Figure 4.3 The fraction of l-fold Ge and l-fold Te atoms (l=1, 2, 3, 4, 5 or 6) as computed from the models simulated with BLYP and PBE.

Figure 4 . 4

 44 Figure 4.4 The inter-polyhedral Ge -Te -Ge (upper panel), the intra-polyhedral Te -Ge -Te (middle panel) and the Te -Te -Te (bottom panel) bond angle distributions. The PBE results are on the left and BLYP results are on the right. Black lines: no account of vdW corrections, red lines: account of vdW corrections.

  [4.5] : Chapter 4. Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4

Figure 4 . 5 4 . 4 .

 4544 Figure 4.5 Top figure: amorphous GeTe 4 total structure factor. Bottom figure amorphous GeTe 4 total pair correlation function. Experimental results from ref. [5] (black dots), BLYP+vdW results (red solid lines) and PBE+vdW results (green solid lines).

Figure 4 . 6

 46 Figure 4.6 Partial pair correlation functions of glassy GeTe 4 . From top to bottom: g GeGe (r), g GeTe (r) and g TeTe (r). PBE+vdW simulations (black lines) and BLYP+vdW simulations (red lines).

Chapter 4 .

 4 Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4 136 gives a coordination numbers of Ge and Te equal to 3.99 and 2.31 respectively. The resulting mean coordination number is then equal to 2.64, in good agreement with the 8-N rule ( n = 2.

Figure 4 . 7

 47 Figure 4.7 The fraction of l-fold Ge and Te atoms (l=1, 2, 3, 4, 5 or 6) as obtained from BLYP+vdW and PBE+vdW models.

137 4 . 4 .

 44 Impact of the XC functional BAD Ge-Te-Ge BAD Te-Ge-Te 0 30 60 90 120 150 180 Angle θ (°) BAD Te-Te-Te

Figure 4 . 8

 48 Figure 4.8 The inter-polyhedral Ge -Te -Ge (upper panel), the intra-polyhedral Te -Ge -Te (middle panel) and the Te -Te -Te (bottom panel) bond angle distributions. PBE+vdW simulations (black lines) and BLYP+vdW simulations (red lines).

Figure 4 . 9

 49 Figure 4.9 The local order parameter distribution for BLYP+vdW system (red lines) and PBE+vdw (black lines). The arrows indicate the positions that correspond to a given defective octahedral geometries. A snapshots from PBE+vdw are also provided.

Figure 4 . 10

 410 Figure 4.10 Partial pair correlation functions of GeTe 4 glass. From top to bottom: g GeGe (r), g GeTe (r) and g TeTe (r). T=300K, Model simulated with BLYP (black lines), model simulated with BLYP+vdW corrections (red lines) and simulation with PBE+vdW (green lines).

Figure 4 . 11

 411 Figure 4.11 The intra-polyhedral Te -Ge -Te (top panel) and the inter-polyhedral Ge -Te -Ge (bottom panel). T=300K, Simulation BLYP@300K (black lines), simulation BLYP+vdW@300K (red lines) and simulation PBE+vdW@300K (green lines).

Figure 4 .

 4 Figure 4.12 The fraction of l-fold Ge and Te atoms (l=1, 2, 3, 4, 5 or 6) as obtained from BLYP@300K (black color), BLYP+vdW@300K (red color) and PBE+vdW@300K (green color) models.

Figure 5 . 3

 53 Figure 5.3 Fragments of the local structure of GST around Ge atoms in the crystalline (left)and amorphous (right) states. Ref[24].

Figure 5 . 5

 55 Figure 5.5 Ab initio molecular-dynamics model of the amorphous state of the phase-change non-volatile memory material Ge 2 Sb 2 Te 5. Ref[28]

Figure 5 . 6

 56 Figure 5.6Crystallization temperatures T x as a function of composition in the Ga-Sb-Te phase phase diagram. Ref[42]

Figure 6 . 1 Figure 6 . 2

 6162 Figure 6.1 Snapshots of crystalline Ga 4 Sb 6 Te 3 model, top figures, from the left to the right;(100), (010) and (001) side views. Bottom figures: 3D snapshots.

Figure 6 . 3

 63 Figure 6.3 Time evolution of the diffusion coefficients of a given chemical species. The insetshows the MSD for Ga, Sb and Te.

Figure 6 . 8

 68 Figure 6.8 Local order parameter distribution for the 3-and 4-fold Ga, Sb or Te atoms. Model A (black lines), model B (red lines) and model C (blue lines).

Figure 6 .Chapter 6 .Figure 6 . 11

 66611 Figure 6.10 (upper plot) Snapshots of C model. Only Sb atoms (yellow) and Te atoms (green) are shown. (bottum plot) Snapshot of the 97-atom Sb-cluster in the 299-atoms model of Ga 4 Sb 6 Te 3 . The atoms of the cluster are shown as black spheres while all the other atoms are represented by gray (red), black and light gray (cyan) lines for Ga, Sb and Te respectively. The shape of the cluster is depicted by the (blue) transparent surfaces of the Voronoi polyhedra built around the Sb atoms of the cluster. Four 299-atoms supercells are shown but the periodic boundary condi-tions are not applied to the atoms of the Sb cluster to highlight its shape and connectivity 177

Figure 6 . 14

 614 Figure 6.14 The proportions of l-fold Ga, Sb and Te atoms (l=1, 2, 3, 4 and 5) as obtained for models D, E and F.

Figure 6 . 15

 615 Figure 6.15 Volume evolution during the quench.

Figure 6 . 16

 616 Figure 6.16 Thermal cycle: the temperature evolution during the simulation.

70 NF= 70 2 N

 70702 (Q) and73 76 F(Q), the Ge -Ge correlations can be eliminated by forming the first difference VII B-Hints on the diffraction theory VIII function ∆F Se (Q) defined as:∆F Se (Q) ≡ 73 76 F(Q) -D(S GeSe (Q) -1) + E(S SeSe (Q) -1) (6.4)where D≡2c Ge c Se (b 73 Ge b 76 Se -Ge b N Se ) and E≡c 2 Se Se ). The Se -Se correlations can be removed by forming the first difference function ∆F Ge (Q) where:

73 73 70 N=

 737370 ) with F≡2c Ge c Se (b 70 Ge b N Se -Ge b 76 Se ) and G≡c 2Ge Ge ). Finally, the Ge-Se correlations can be removed by forming the first order difference function ∆F Se ′ (Q) defined as:∆F Se ′ (Q) ≡ 73 76 F(Q) -b 73Ge b 76 Se b 70 Ge b N Se H(S GeSe (Q) -1) + I(S SeSe (Q) -1) (6.8) where H≡c 2 Ge (b 2 73 Ge -b 73 Ge b 76 Se b 70 Ge b N Se b 2 70 Ge ) and I≡c 2 Se (b 2 76 Se -b 73 Ge b 76 Se b 70 Ge b N Se b 2 N Se ).

  = 10.0(1), b 73 Ge = 5.09(9) and b 76 Se = 12.2(1) fm [2]. The atomic fraction of Ge and Se are c Ge = 1/3 and c Se = 2/3 respectively.

on the diffraction theory 5 -( 6 . 29 ) 6 - 1 -

 562961 From Bhatia-Thornton to Ashcroft-LangrethS 11 (Q) = c 1 S NN (Q) + S CC (Q)/c 1 + 2S NC (Q) (6.26) S 22 (Q) = c 2 S NN (Q) + S CC (Q)/c 2 -2S NC (Q) (6.27) S 12 (Q) = (c 1 c 2 ) 1/2 S NN (Q) -S CC (Q)/(c 1 c 2 ) 1/2 + (c 2 /c 1 ) 1/2 -(c 1 /c 2 ) 1/2 S NC (Q) (6.28) From Bhatia-Thornton to Faber-Zimana 11 (Q) = S NN (Q) + c 2 1 S CC (Q) + (2/c 1 )S NC (Q) -c 2 /c 1 (6.30) a 22 (Q) = S NN (Q) + c 2 2 S CC (Q) -(2/c 2 )S NC (Q) -c 1 /c 2 (6.31) a 12 (Q) = S NN (Q) -S CC (Q)/(c 1 c 2 ) + ( Paircorrelation functionThe pair correlation function (PCF) allows access to the real space atomic arrangements. This quantity expresses the probability of finding a neighbor in a shell dr at the distance r of a given atom, Figure.[6.17].

Figure 6 . 17

 617 Figure 6.17 Real space discretization for g(r) calculation.

Figure 6 . 18

 618 Figure 6.18 Coordination numbers in the case of a GeSe 2 system.

Figure 6 . 19

 619 Figure 6.19 Illustration of the shortest-path criterion for ring statistics: in this specific example, starting from a node, three possible rings could be found. For each ring the bond between the node and the two neighbors (from which the research began) are drawn in solid line, the rest of the path is drawn in dashed lines.

Figure [ 6 . 20 ]Figure 6 . 20

 620620 Figure 6.20 Left figure: Two-dimensional illustration of the Voronoi regions in a system of balls of different sizes. The edges of the Voronoi polyhedra are shown by thick lines. The dashed lines depict the edges of the conventional Voronoi polyhedra of the system as computed for the centers of the balls. Right figure: Different types of configurations of four balls in 3D (left) and three disks in 2D (right). (a) Singlet set: one inscribed sphere; (b) Double set: two inscribed spheres; (c) Nonsimpliciable: no inscribed spheres. This figure contains two figures shown as separate in ref. [5].

Figure 6 . 21

 621 Figure 6.21 Snapshots of tetrahedral, octahedral and defective octahedral configurations.The value of the local order parameter is provided.

Figure[ 6 .Figure 6 .

 66 Figure[6.22] is a snapshot of a three-dimensional view of structural subunits of amorphous GeSe 2 highlighting different type of connections:

  the original ionic positions are expressed as a function of h as: R I = hS I (6.53) XXI E-Parrinello-Rahman barostat XXII

  Ω φ * i (r)φ i (r)dr = Ω Ω unit φ * i (hs)φ i (hs)ds = Ω unit φ * i (s)φ i (s)ds(6.55)

  [1.43]) written for the scaled wavefunctions and with and additional constant term pΩ. The corresponding equations of motion us -pδ us (h T )

  n(r) = ∑ µν P µν ϕ µ (r)ϕ ν (r) (6.63) XXV F-GPW basis set XXVI with P µν is a density matrix element, ϕ ν (r) = ∑ i d iµ g i (r) with g i (r) are a primitive Gaussian functions and d iµ is their corresponding contraction coefficients.The auxiliary plane wave basis set ( ñ(r)) is introduced as:n(r) = ñ(r) + ∑ A n A (r) -∑ A ñA (r)(6.64)ñ(r) is the soft contribution to the density from the interaction regions, expanded in a plane wave basis set as follows: ñ(r) = 1 Ω ∑ G ñ(G)exp(iG.r) (6.65)

Figure 6 . 23 2 ∇ 2

 62322 Figure[6.23] shows an illustration of the wave function expansion with the GPW method. Near the the atom A the density is given by the density of the atom A expanded in the localised Gaussian basis set (n(r) = n A (r)). However, far from the atom the density is equal to the one expanded on a plane wave basis set. The way the density is expanded in the GPW method allow the efficient and rapid conversion between n(r), ñ(r) and ñ(G) by mean of FFT and a reliable mapping procedure. Within this representation the KS DFT energy (Eq.1.13 in section [1.2])

Figure 6 . 24

 624 Figure 6.24 Partial pair correlation functions as computed with HSE06 and BLYP functionals.

❅❅

  Phd related papers: ❅ Note: Accounting for pressure effects on the calculated equilibrium structure of glassy GeSe 2 . J. Chem. Phys. 137, 046101 (2012); Assil Bouzid and Carlo Massobrio. ❅ The density-driven defect-mediated network collapse of GeSe 2 glass. Phy. Rev. B 90, 054206 (2014); Kamil Wezka, Assil Bouzid, Keiron J. Pizzey, Philip S. Salmon, Anita Zeidler, Stefan Klotz, Henry E. Fischer, Craig L. Bull, Matthew G. Tucker, Mauro Boero, Sébastien Le Roux, Christine Tugéne and Carlo Massobrio. ❅ First principles study of the amorphous Ga 4 Sb 6 Te 3 phase change compound. In preparation; Assil Bouzid, Silvia Gabardi, Mauro Boero, Carlo Massobrio and Marco Bernasconi. Performances of the XC functionals and effects of van der waals forces on the structure of glassy GeTe 4 . In preparation. ❅ First principles molecular dynamics study of GeSe 4 glasse under pressure. In prepa-Structural properties of glassy Ge 2 Se 3 from first-principles molecular dynamics. Physical review B 86, 224201 (2012); Sébastien Le Roux, Assil Bouzid, Mauro Boero, and Carlo Massobrio.

GeS 2 :

 2 Atomic structure and bonding properties. Physical review B 88, 174201 (2013); Massimo Celino, Sébastien Le Roux, Guido Ori, Benoit Coasne, Assil Bouzid, Mauro Boero and Carlo Massobrio.

  GeS 2 : A realistic model based on a first-principles approach. Phys. Rev. B 90, 045423 (2014); Guido Ori, Carlo Massobrio, Assil Bouzid, Mauro Boero, Massimo Celino, Benoit Coasne.

  xxx 10 Configurations atomistiques du verre Ga 4 Sb 6 Te 3 avec les principaux motifs structuraux. . . . . . xxxi 1.1 General scheme of a computer code for MD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.2 Schematic diagram of "Jacob's ladder" of exchange correlation functionals proposed by J. P. Perdew. 20 1.3 Schematic illustration of all-electron (dashed lines) and pseudopotential (solid lines) and their corresponding wave functions. The radius at which the all-electron and pseudo-electron values match is designated as r c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.

1

Fonctions de corrélation de paires partielles pour l'amorphe GeSe 2 . Voir

[3.1]

pour plus de détails. xxi 2 Cycle thermique: l'évolution de la température au cours de la simulation. . . . . . . . . . . . . . xxii 3 La variation en fonction de la pression des fractions d'atomes de Ge "Corner-Sharing" et "Edgesharing", ainsi que la répartition de ce dernier dans des configurations de type Ge 1 , Ge 2 and Ge 3&4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii 4 Configurations atomistiques du verre GeSe 2 prises à différentes pressions. Couleur violette: Atomes de Ge et couleur jaune: atomes de Se. . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv 5 Fonctions de corrélation de paires partielles pour l'amorphe GeSe 4 sous pressions. Du haut jusqu'en bas: g SeSe (r), g GeSe (r), g GeGe (r) et fonction de corrélation radiale totale g tot (r). . . . . . xxvi 6 Configurations atomistiques du verre GeSe 4 à pression ambiante et à 13 GPa. La flèche bleue indique une configuration de type Se 3 , la rouge une configuration de type Ge 5 et la verte indique un défaut de type Ge -Ge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxvii 7 Fonctions de corrélation de paires partielles pour l'amorphe GeTe 4 . Du haut jusqu'en bas: g GeGe (r), g GeTe (r) et g TeTe (r) De gauche à droite: modèle généré avec PBE et modèle généré avec BLYP. Simulations avec les corrections de vdW (lignes rouges) et les simulations sans les corrections de vdW (lignes noires). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxviii 8 La distribution de paramètre d'ordre local pour le modéle BLYP+vdW (couleur rouge) et pour le modéle PBE+vdw (couleur noire). Les flèches indiquent les positions qui correspondent à des configurations octaédriques défectueuse. Des photos du modéle PBE+vdw sont également fournis. xxix 9 L'évolution de la température au cours de la trempe du modèle A (couleur verte en pointillés), du modèle B (couleur bleue) et du modèle C (couleur rouge). . . . . . . . . . . . . . . . . . . . . . .

  2 at ambient pressure for various convergence criteria and various time steps. Top panel: red line: CP, 5 a.u.; black dashed line: CP, 10 a.u. Middle panel: black dashed line: CP, 10 a.u.; green line: BO, 10 a.u., 10 -6 a.u.; magenta line: BO, 100 a.u., 10 -6 a.u.. Bottom panel: black dashed line: CP, 10 a.u.; magenta line: BO, 100 a.u., 10 -6 a.u.; blue line: BO, 100 a.u., 10 -5 a.u.. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Car-Parrinello conserved total energy variation for different temperatures. Tests performed on a GeSe 2 system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 The simulation time per molecular dynamic step variation as a function of number of CPUs used.

(a): a system of Ga 4 Sb 6 Te 3 with 117 atoms has been simulated on IBM x3750 M4 (Ada-IDRIS). (b) results of Ga 4 Sb 6 Te 3 simulation with 299 atoms on SGI ICE machine (Jade-CINES). (c) results of Ga 4 Sb 6 Te 3 simulation with 299 atom on IBM Blue Gene / Q (Turing-IDRIS) . . . . . . . . . . 1.8 Temperature variation and electronic fictitious kinetic energy as imposed by the velocity rescaling method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Schematic representation of

  . . . . . . . . 2.10 Representative subset of Ge and Se atoms in amorphous GeSe 4 where Ge atoms are dark (blue) and Se atoms are light (green). Se atoms along a connection path between two Ge atoms are labelled as AA, those between one Ge atom and one Se atom are labelled as AB, and those between two Se atoms are labelled as BB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 The calculated partial pair correlation functions for glassy GeSe 4 . The label Se-Se indicates g

SeSe (r), the label Ge-Se indicates g GeSe (r) and the label Ge-Ge indicates g GeGe (r). The circle symbols curve corresponds to PW, the other one to PBE. A set of representative error bars are given in the g GeGe (r) case. Figure taken from Ref.

[32]

. . . . . . . . . . . . . . . . . . . . . . . 2.12 Pressure-volume equation of state for amorphous GeSe 2 . Experimental compression in a hydrostatic medium, a 4:1 methanol:ethanol mixture (solid circles); experimental decompression from 8.55 GPa (open triangles); simulation data open squares Ref.

  .10 The pressure dependence of the fractions of Ge CS atoms and ES Ge atoms is given, together with a breakdown of the latter into its contributions from Ge 1 , Ge 2 and Ge 3&4 units. Left plot: NVT simulation results, right plot: NPT simulation results. . . . . . . . . . . . . . . . . . . . . . . . .

[(magenta)]; or (iii) FPMD NVT simulations [broken (black) curve], or cold compression NVT simulations at T= 300 K [dotted broken (maroon) curve with cross marks], or NPT simulations [dotted (dark green) curve with stars] or HSE06 simulation [orange cross mark]. In the bottom figure, the horizontal chained line gives the

"8-N" rule expectation of n = 2.67. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

  plot: NVT simulation results, right plot: NPT simulation results. . . . . . . . . . . . . . 3.16 Atomistic configurations taken from GeSe 2 glass at different pressures. Ge atoms are dark (purple) and Se atoms are light (yellow). Bonds are drawn when two atoms are separated by a cutoff distance given by the position of the first minimum in g GeSe (r). . . . . . . . . . . . . . . . . . . . 3.17 Pressure-Volume equation of state for amorphous GeSe 4 . Solid black circles are experimental data from ref.

  .20 Pair correlation function of amorphous GeSe 4 under pressures. from top to bottom g SeSe (r),

(black lines with diamond). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 The pressure evolution of the partial Faber-Ziman structure factor from from Fourier transform of the partial real space counterpart presented on figure.[3.20]. From the top to the bottom S FZ SeSe (k), S FZ GeSe (k) and S FZ SeSe (k). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

  (black dots), BLYP results (blue solid lines), BLYP+vdW results (red solid lines), PBE results (purple dashed lines) and PBE+vdW results (green solid lines). PBE and PBE+vdW was shifted to y+2 (for the s(k)) and to y+2.5 (for the g(r)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Partial pair correlation functions of GeTe 4 glass. From top to bottom: g GeGe (r), g GeTe (r) and g The fraction of l-fold Ge and l-fold Te atoms (l=1, 2, 3, 4, 5 or 6) as computed from the models simulated with BLYP and PBE. . The local order parameter distribution for BLYP+vdW system (red lines) and PBE+vdw (black lines). The arrows indicate the positions that correspond to a given defective octahedral geometries. A snapshots from PBE+vdw are also provided. . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Partial pair correlation functions of GeTe 4 glass. From top to bottom: g GeGe (r), g GeTe (r) and g

TeTe (r). From left to right: PBE model and BLYP model. Simulations with vdW corrections (red lines) and simulations without vdW corrections (black lines). . . . . . . . . . . . . . . . . . . . . 4.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 The inter-polyhedral Ge-Te-Ge (upper panel), the intra-polyhedral Te-Ge-Te (middle panel) and the Te -Te -Te (bottom panel) bond angle distributions. The PBE results are on the left and BLYP results are on the right. Black lines: no account of vdW corrections, red lines: account of vdW corrections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Top figure: amorphous GeTe 4 total structure factor. Bottom figure amorphous GeTe 4 total pair correlation function. Experimental results from ref. [5] (black dots), BLYP+vdW results (red solid lines) and PBE+vdW results (green solid lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Partial pair correlation functions of glassy GeTe 4 . From top to bottom: g GeGe (r), g GeTe (r) and g TeTe (r). PBE+vdW simulations (black lines) and BLYP+vdW simulations (red lines). . . . . . . . 4.7 The fraction of l-fold Ge and Te atoms (l=1, 2, 3, 4, 5 or 6) as obtained from BLYP+vdW and PBE+vdW models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 The inter-polyhedral Ge-Te-Ge (upper panel), the intra-polyhedral Te-Ge-Te (middle panel) and the Te -Te -Te (bottom panel) bond angle distributions. PBE+vdW simulations (black lines) and BLYP+vdW simulations (red lines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 TeTe (r). T=300K, Model simulated with BLYP (black lines), model simulated with BLYP+vdW corrections (red lines) and simulation with PBE+vdW (green lines). . . . . . . . . . . . . . . . .

  Schematic presentation of the possible crystallization process in amorphous Ge 2 Sb 2 Te 5 and GeTe.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4

  Volume evolution during the quench. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV 6.16 Thermal cycle: the temperature evolution during the simulation. . . . . . . . . . . . . . . . . . . V 6.17 Real space discretization for g(r) calculation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . XIII 6.18 Coordination numbers in the case of a GeSe 2 system. . . . . . . . . . . . . . . . . . . . . . . . . XIV 6.19 Illustration of the shortest-path criterion for ring statistics: in this specific example, starting from

a node, three possible rings could be found. For each ring the bond between the node and the two neighbors (from which the research began) are drawn in solid line, the rest of the path is drawn in dashed lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XVI 6.20 Left figure:

  Percentages of Te and Ge atoms belonging to Te and Ge bonds, N TeTe and N GeGe . First maxima and minima of the partial pair correlation functions and the partial coordination number n αβ in the amorphous and liquid. Also, coordination numbers of different species, mean coordination number in the amorphous and liquid and proportion of homopolar bonds for different species. For g TeTe (r) the cutoff distance for computing the coordination numbers was fixed to 3.2 Å.6.2 First maxima and minima of the partial PCF and the partial coordination number n αβ for Models A, B and C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Inter-atomic distances as found in GaTe Crystal. Experimental results are compared to CPMD-GTH and CPMD-MT results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Partials and coordination numbers for models A, B and C. . . . . . . . . . . . . . . . . . . . . . 6.5 Nearest-neighbor analysis of model A, B and C of amorphous Ga 4 Sb 6 Te Partials and mean coordination numbers for models D, E and F. . . . . . . . . . . . . . . . . . . 6.8 The weighting coefficients (barns (10 -28 m 2

	Also shown the
	percentages of Ge atoms found in edge-sharing configurations (N ES Ge ) and corner-sharing configu-rations (N CS . .
	4.3 r i j : position of the first maximum of the pair distribution functions g i j ), partial coordination num-
	bers n i j (integration of g i j up to the first minimum), coordination numbers n

. . . . . . . . . . . . . . . 3.3 NVT pressure values at which comparisons of the simulation and experiment data sets are made. NPT simulations data are also provided. The corresponding number densities and reduced densities are given. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Pressure points at which a comparison between simulation and experiments is performed. . . . . . 4.1 Bond distances r i j (position of the first maximum of the pair distribution functions g i j ), partial coordination numbers n i j (integration of g i j up to the first minimum ), coordination numbers n i and the mean coordination number n. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Ge ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i and the mean coordination number n. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Chronology of phase-change materials developments [2, 20]. . . . . . . . . . . . . . . . . . . . . xiii 6.1 3 , only configurations that contribute to more than 2% are shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 First maxima and minima of the partial PDFs and the partial coordination numbers n αβ in the amorphous phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7

  1-x , une autre composition a été le sujet d'une étude sous pression. Le verre GeSe 4 correspond à la concentration x=0.2 qui existe à la limite de la fenêtre des concentrations où les matériaux présentent une certaine flexibilité, et ceci à cause de la présence de nombreuses chaines de Se. Les simulations ont été effectuées dans l'ensemble NPT, permettant ainsi de fixer une valeur cible au début du calcul. Pour couvrir un domaine de pressions allant de l'ambiante jusqu'à 14 GPa, 5 systèmes indépendants ont été générés et chaque système à été l'objet d'un cycle thermique comme dans le cas du GeSe 2 .

  13.2 GPa.

		13.2 GPa
	(r) SeSe	8.18 GPa 6.7 GPa
	g	4.16 GPa
		0.15 GPa
	GeSe (r)	6.7 GPa 8.18 GPa 13.2 GPa
	g	4.16 GPa
		0.15 GPa
	GeGe (r)	6.7 GPa 8.18 GPa 13.2 GPa
	g	4.16 GPa
		0.15 GPa
		13.2 GPa
	g tot (r)	6.7 GPa 8.18 GPa
		4.16 GPa
		0.15 GPa
		2	4	6	8	10
			Distance r (Å)
						xxv

  deux états garantira aussi une conservation des données sur des longues périodes. Aussi, un large contraste dans les propriétés optiques et/ou de résistivité permettra la lecture du bit d'information enregistré. En dernier lieu, pour pouvoir écrire et effacer nos données, le PCM utilisé doit avoir un nombre aussi grand que possible de cycles réversibles.Récemment, le Ga 4 Sb 6 Te 3 (GaST) a été proposé comme une alternative possible du PCM couramment utilisé, le Ge 2 Sb 2 Te 5

Système ternaire: Ga 4 Sb 6 Te 3 La deuxième partie de la thèse est dédiée à l'étude des matériaux à changement de phase (PCM). Ces matériaux se caractérisent par la transition réversible et rapide entre l'état cristallin et l'état amorphe. Ceci permet l'utilisation des PCM dans la fabrication des dispositifs de stockage de données. Trouver un bon PCM qui satisfait les propriétés optimales de la mémoire universelle est loin d'être simple. Un bon PCM doit avoir une grande vitesse de transition entre l'état amorphe et cristallin permettant ainsi une grande vitesse d'écriture de données. Une xxix grande stabilité des

Table 1 .

 1 1.6. CPMD implementation ambient conditions. The following table illustrates the different tests results. The outcome of these tests over a significant time scale of 10ps are displayed on figure[1.5],

	Run Method Time step (a.u.) Convergence (a.u.) Conservation (a.u./ps) Time (h)
	CP1 CP2 BO1 BO2 BO3	CP CP BO BO BO	5 10 10 100 100	--10 -6 10 -5 10 -6	5.2 × 10 -9 1.1 × 10 -8 2.6 × 10 -4 2.1 × 10 -4 6.1 × 10 -6	2.27 1.05 37.58 4.1 6.94

1 Energy conservation in a.u./ps for Car-Parrinello (CPMD) and Born-Oppenheimer (BOMD) molecular dynamics tests and CPU timing tested for a system of GeSe 2 at ambient pressure. Results are displayed for a 10 ps of dynamics on LINUX cluster (Intel processors Xeon L5420) using 48 processors for all tests, the runs were performed using the 3.17 CPMD version.

  -5 a.u. (run BO2). Unfortunately (for a time step 100 a.u.), the results show a large deviation from the stable CP1 energy, with an energy conservation of 2.1 × 10 -4 a.u./ps.

		1e-06							(a)
		0
		-2e-06 -1e-06						CP1 CP2
	(a.u)	0e+00 5e-04							(b)
	E cons	-1e-03 -5e-04						CP2 BO1 BO3
		2e-03 4e-03						CP2 BO3 BO2	(c)
		0e+00
		-2e-03	0	2	4	Time (ps)	6	8	10

Figure

1

.5 Conserved total energy E cons as computed from Eq.

[1.46]

from Car-Parrinello (CPMD) and Born-Oppenheimer (BOMD) molecular dynamics for a system of GeSe 2 at ambient pressure for various convergence criteria and various time steps. Top panel: red line: CP, 5 a.u.; black dashed line: CP, 10 a.u. Middle panel: black dashed line: CP, 10 a.u.; green line: BO, 10 a.u., 10 -6 a.u.; magenta line: BO, 100 a.u., 10 -6 a.u.. Bottom panel: black dashed line: CP, 10 a.u.; magenta line: BO, 100 a.u., 10 -6 a.u.; blue line: BO, 100 a.u., 10 -5 a.u..

  illustrate the outcome of the tests:

	0,5 1 2 4 1 2 4 8 16 32 64 2 4 8 16 32 64 128 CPU time / MD step CPU time / MD step CPU time / MD step	128 512 64	128	256	512 Number of CPUs 1024 256 512	[MPI x 4 OMP] 1024 [MPI x 2 OMP] 2048 1024 [MPI x 2 OMP] [MPI x 4 OMP] [MPI x 8 OMP]

Table 3 .

 3

	Simulation at ρ exp	1]	20 27	30 19	58 54	22 19
	Simulation at ρ b		23	18	35	42
	Simulation Ref. [6]		17	30	38	45
	Experiment: Ref. [5]		25	20	34	41

2 N Ge-Ge (N Se-Se ) is the percentage of Ge (Se) atoms in Ge-Ge (Se-Se) homopolar bonds, N Ge (ES) is the percentage of Ge atoms forming edge-sharing connections and N Ge (CS) is the percentage of Ge atoms forming corner-sharing connections.

Table 3 .

 3 The pressure-volume equation of state for GeSe 2 glass under compression where

		03262	1	0	0.0334	1	0	0.0324	1
	0.37	0.0334	1.0239	-	-	-	-	-	-
	1.29	0.03505	1.0745	-	-	-	-	-	-
	2.15	0.03711	1.1376	-	-	-	-	-	-
	3.55	0.03929	1.2045	3.0(5)	0.0396	1.1856	-	-	-
	4.88	0.04123	1.2639	4.7(5)	0.0428	1.2814	-	-	-
	6.87	0.04395	1.3473	-	-	-	-	-	-
	7.45	0.04453	1.3651	6.3	0.0455	1.3622	-	-	-
	8.73	0.04639	1.4221	7.1(5)	0.0469	1.4041	-	-	-
	9.87	0.04771	1.4626	8.2(5)	0.0489	1.4640	8.7	0.0475	1.4660
	11.56	0.04985	1.5282	-	-	-	10.9	0.0499	1.5401
	13.6	0.05139	1.5754	-	-	-	12.8	0.0519	1.6018
	15.27	0.05347	1.6392	-	-	-	14.4	0.0534	1.6481
	-	-	-	-	-	-	16.1	0.0549	1.6944
				NPT Simulation
		P(GPa)			ρ sim (Å -3 )			ρ sim /ρ 0
		0.13			0.03225			1
		4.05			0.03825			1.1858
		5.46			0.04003			1.2410
		9.48			0.04515			1.4000
		16.21			0.05277			1.6362

3 NVT pressure values at which comparisons of the simulation and experiment data sets are made. NPT simulations data are also provided. The corresponding number densities and reduced densities are given. P denotes the pressure, V is the volume at pressure P and V 0 is the volume under ambient conditions. The measured data from Mei and Co. ((black) • with vertical error bars) are compared to the results obtained from first-principles molecular dynamics in the NVT present work ((red) ), NPT present work results ((magenta) ×) and in the work by Durandurdu and Drabold ((blue) ♦). The computed and the measured data are fitted to a second-order Birch-Murnaghan equation of state (solid blue and solid black curves respectively).

Table 3 . 4

 34 Pressure points at which a comparison between simulation and experiments is performed.

	:

Figure 3.17 Pressure-Volume equation of state for amorphous GeSe 4 . Solid black circles are

experimental data from ref.

[23]

, red square are simulation data. Green-dashed line represent the third-order Birch-Murnaghan fit to the experimental data.

  result from Fourier transforming the real space partial pair distribution functions. At ambient pressure one notices that the FSDP stems

	115				3.2. Amorphous GeSe 4 under pressure
	2,5							4
	7 1 1,5 2 1 S ) Peak position (Å -1 (k) SeSe FZ	(b) (a) 0 0	2 2	4 4	6 Pressure (GPa) 8 6 8 Exp-PP height Sim-PP height	10 10	12 Sim-FSDPP Exp-PPP Exp-FSDPP Sim-PPP fromd direct 14 Sim-PPP from FFT 12 14 0 1 2 3 13.2 GPa PP height 8.18 GPa 6.7 GPa 4.16 GPa 0.15 GPa
	6							13.2 GPa
	3 4 5 1 GeSe S X (Q) S FZ (k)						6.7 GPa 8.18 GPa 13.2 GPa	4.16 GPa 0.15 GPa 6.7 GPa 8.18 GPa
	2						4.16 GPa
	0 1 0 GeGe (k) S FZ	2	4 Scattering wave vectors k (Å 6 8 10 12 -1 ) 0 GPa 14	16 0.15 GPa 4.16 GPa 6.7 GPa 8.18 GPa 13.2 GPa
	Figure 3.18 (a) The FSDP position (FSDPP) as obtained from direct calculation of the total structure factor (green triangle up) and experiments (solid black circles). The principle peak position (PPP) as computed from direct calculation of the total structure factor (magenta triangles down), from Fourier transformed total struc-ture factor (green left triangles), and experiments (blue squares). In addition the 2 4 6 8 10 12 14 Scattering wave vectors k (Å -1 )
	PPP height are provided: as computed from Fourier transformed total structure
	factor (red plus) and experiment (black stars). (b) The pressure evolution of the

X-ray structure factor from FPMD direct calculation on reciprocal space (light green lines), from Fourier transform of the partial real space counterpart (blue dashed lines) and X-ray diffraction from Ref.

[23]

(black lines with diamond).

table : Chapter 4 .

 :4 Performances of the XC functionals and effects of Van der Waals forces on the structure of GeTe 4 130 r GeGe (Å) r GeTe (Å) r TeTe (Å) n GeGe n GeTe n TeTe n Ge n Te

					n
	PBE	2.67	2.65	2.88	0.26 1.06 2.71 4.53 3.78 3.93
	PBE+vdW	2.49	2.65	2.83	0.15 1.04 2.46 4.34 3.51 3.68
	BLYP	2.46	2.62	2.89	0.28 0.96 1.60 4.14 2.56 3.14
	BLYP+vdW	2.43	2.59	2.84	0.37 0.90 1.40 3.99 2.31 2.64

Table 4 . 1

 41

Table 4 . 2

 42 Percentages of Te and Ge atoms belonging to Te and Ge bonds, N TeTe and N GeGe .

	Ge

  r TeTe (Å) n GeGe n GeTe n TeTe n Ge n Te n

	Blyp@300K	2.49	2.62	2.78	0.20 0.96 2.35 4.06 3.31 3.46
	Blyp+vdW@300K	2.47	2.60	2.78	0.16 0.96 2.03 4.02 3.0 3.20
	PBE+vdW@300K	2.54	2.62	2.83	0.21 1.00 2.49 4.23 3.49 3.64
	Table

Table 6 . 1

 61 First maxima and minima of the partial pair correlation functions and the partial coordination number n αβ in the amorphous and liquid. Also, coordination numbers of different species, mean coordination number in the amorphous and liquid and proportion of homopolar bonds for different species. For g TeTe (r) the cutoff distance for computing the coordination numbers was fixed to 3.2 Å.

	displays the different pair correlation

Table 6 . 2

 62 First maxima and minima of the partial PCF and the partial coordination number n αβ for Models A, B and C.

	GPW, see appendix D for

Table 6 . 3

 63 Inter-atomic distances as found in GaTe Crystal. Experimental results are compared to CPMD-GTH and CPMD-MT results.

Table 6 . 4

 64 . Partials and coordination numbers for models A, B and C.

		model A model B model C
	n Ga n Sb n Te n	3.31 3.18 2.55 3.08	3.79 3.28 2.80 3.33	3.96 3.31 2.80 3.39

Table 6 . 6

 66 First maxima and minima of the partial PDFs and the partial coordination numbers n αβ in the amorphous phase.

	table:

Table 6 . 7

 67 The coordination numbers are given in the following table: Partials and mean coordination numbers for models D, E and F.Interestingly, model D features similar coordination numbers as Model B and Model C. Ga is almost 4-fold, Sb 3-fold and Te is in a mixture of 2-and 3-fold configurations. For more details, the fraction of l-fold Ga, Sb and Te atoms for all quenched models is given below:

		model D model E model F
	n Ga n Sb n Te n	3.79 3.29 2.80 3.33	3.87 3.33 2.76 3.36	3.84 3.27 2.69 3.31

Table 6 . 8

 68 The weighting coefficients (barns (10 -28 m 2 )) for Ge-Ge, Ge-Se and Se-Se Faber-Ziman partial structure factors in the total and difference functions (eq.[6.2] to eq.[6.7]). The numerical values take into account the isotopic enrichments of the samples used in the experiments, the coherent scattering length are given by: b 70 Ge

	70 N F(Q)	0.1111(22)	0.354(4)	0.2823(6)
	N N F(Q)	0.0744(4)	0.2899(8)	0.2823(6)
	73 76 F(Q)	0.0288(5)	0.276(3)	0.662(11)
	∆F Se (Q)	0.0000(7)	0.184(8)	0.588(13)
	∆F Ge (Q)	0.099(7)	0.237(9)	0.0000(13)
	∆F Se ′ (Q)	-0.057(18)	0.0000(18)	0.442(18)

  2n SeGe(6.37) n Ge = n GeGe + n GeSe(6.38)nSe = n SeSe + n SeGe (6.39) n = c Ge * n Ge + c Se * n Se (6.40) (6.41) In the case of ternary system: say, Ga 4 Sb 6 Te 3 , we have: Accordingly, the mean cordination number is calculated as follws: n = c Ga * n Ga + c Sb * n Sb + c Te * n Te

	C-Analysis techniques and definitions	XVI
	4n GaSb = 6n SbGa	(6.42)
	4n GaTe = 3n TeGa	(6.43)
	6n SbTe = 3n TeSb	(6.44)
	n Ga = n GaGa + n GaSb + n GaTe	(6.45)
	n Sb = n SbGa + n SbSb + n SbTe	(6.46)
	n Te = n TeGa + n TeSb + n TeTe	(6.47)
		(6.48)
		(6.49)
		XV

Acknowledgments

This dissertation is the result of three years of research that I carried out within the atomic scale modelling team in the department of inorganic chemistry, at the IPCMS. Certainly, the work described in this thesis could not have been accomplished without the help of many people. Now, at the end of this adventure it is the perfect time to give my sincere acknowledgements.

Part II

Binary chalcogenides

CHAPTER 2

STATE OF THE ART AND MOTIVATIONS Summary

In this chapter the main steps toward an accurate description of the Ge x Se 1-x chalcogenides will be presented. The main focus will be on the theoretical achievements obtained within the framework of first-principles molecular dynamics. When needed, a close connexion with the experimental findings will be established. Mostly, I shall restrict myself to accomplishment obtained within the IPCMS team, the purpose being to set the scene for the motivation and the layout of my own work.

Introduction

The structure of the binary chalcogenides has attracted a high scientific interest for a long time. In order to describe and understand the atomic scale picture of chalcogenides, many experiments, theories and models have been developed. In particular, atomic-scale simulations are required in order to elucidate the structural properties, a necessary prerequisite to any further study linking the structure to the macroscopic, measurable behavior. On the molecular dynamics side, one can classify the researches performed by the Vashishta team in the late 80's as the first milestone of atomic-scale modelling for chalcogenides. The simulation scheme that Vashishta used was based on classical molecular dynamics employing effective potentials. The outcome of these classical simulations was discussed in section [1.1.2.2]. In what follows, we shall refer almost exclusively to first-principles molecular dynamics simulations fully based on (or at least highly inspired by) density functional theory. On the side of CHAPTER 3

BINARY CHALCOGENIDES UNDER

PRESSURE Summary

In this chapter, results on GeSe 2 and GeSe 4 simulations under pressure will be presented. The real and reciprocal space properties will be discussed and compared to neutron diffraction experiments. Details on the structural transitional mechanisms under pressure will be provided.

Amorphous germanium diselenide under pressure

The initial GeSe 2 configuration was taken from the system studied in ref. [1]. Simulations were performed within the FPMD framework, namely the Car-Parrinello method as implemented in the CPMD code. We adopted a generalized gradient approximation (GGA) for the exchange and correlation part of the total energy. In particular, the Becke-Lee-Yang-Parr (BLYP) has been used. The Troullier-Martins norm-conserving pseudopotentials are employed for the corevalence interactions description. The largest cutoff used in the pseudopotentials construction was equal to 1.06 Å. A fictitious electron mass of 1000 a.u. and a time step of ∆t = 0.24 fs were adopted to integrate the equations of motion, ensuring good control of the conserved quantity. The wave functions are expanded in plane waves basis set at the Γ point of the super cell, with a cutoff energy of 20 Ry. A temperature control is implemented for both ionic and electronic degrees of freedom by using the Nosé-Hoover thermostats. The methods used to generated the arrangement in the pseudo-cubic fashion, we computed the local order parameter distribution. The local order parameter (q-parameter) is defined by:

where θ i jk is the angle formed between a central Ge atom j and its neighboring atoms i and k.

The q-parameter is equal to 0 for a perfectly Ge centered octahedral Ge and 1 for a perfectly Ge centered tetrahedron. Values in between 0 and 1 correspond to the occurrence of defective octahedral configuration (see annex C for more details). Figure. [3.13] illustrate the evolution of q-parameter distribution when the density increases.

Figure 3.13

The reduced density dependence of the q-parameter distribution as obtained from simulations of GeSe 2 glass. At each density, the distribution for all Ge atoms is broken down into its contributions from n-fold coordinated Ge atoms (n = 3, 4, 5 or 6). The vertical (red) arrows mark the q value expected for trigonal bipyramidal units.

At ambient (ρ/ρ 0 = 1) the distribution shows a maximum at q ∼ 0.0966, value close to the crystalline phase one (q ∼ 0.0979) [15], reflecting that most of the Ge atoms are in 4-fold coordinated tetrahedral units. As the pressure increase (ρ/ρ 0 = 1.463), the intensity of this peak reduces and a sharp peak at small negative q-value rises corresponding to the formation of distorted octahedral environments. The breakdown of the q-distribution into the contribution CHAPTER 4

PERFORMANCES OF THE XC FUNCTIONALS AND EFFECTS OF VAN DER WAALS FORCES ON THE STRUCTURE OF

GeTe 4

Summary

In this chapter, the structure of glassy GeTe 4 is studied, in particular the effects of the account of the dispersion forces on the structure of this compound is detailed. In addition, comparison of the performances of BLYP and PBE XC functional is presented.

Introduction

Van der Waals (vdW) forces result from the sum of the attractive or repulsive forces between molecules other than those due to Coulomb interactions, hydrogen bonds, or electrostatic interaction between ions. Bonds due to vdW interactions differ from ionic and covalent bonding. They are generated by correlations between the fluctuating polarizations of nearby particles.

Even though bonds due to vdW interactions are very weak compared to covalent bonds, they play an important role in the cohesive energy of all kinds of materials: intermolecular complexes, extended systems, liquids, or biological macromolecules. As a consequence, the account of these forces on any computational scheme is highly recommended. In this thesis, we have

STATE OF THE ART Summary

In this chapter we recall the essence on phase change materials and the physics behind their outstanding behavior. In particular, we give insight into the performances that a universal memory should feature in addition to the reasons behind the high crystallisation speed of PCMs. Previous experimental works and motivations of the present work are presented.

Introduction

Since the dawn of mankind, humans tried to save knowledge and chronicle events by painting and inventing scripture. The very early one being the hieroglyphs, invented by the Egyptians. Later on, as the civilization progresses, men discovered and invented more and more efficient and sophisticated storage tools. The 20 th century was marked by the plethora of new advanced techniques able to save huge amount of data. In particular, the non-volatile 'flash' memory (flash-RAM) [1] is one of the most used memory device, the technology behind relying on the magneto-optical properties of the materials employed. During the past decades, three generations of optical storage products were developed and optimised, the compact disks (CD), digital versatile disks (DVD) and Blue-ray disks (BD). These memory devices find applications in the microelectronics industry, for example in the smart-phones, personal-computers and digital cameras.

Despite all these progresses, the endless advances in new technologies such as the forth generation networks, the automobile new concepts, the data mining etc require new materials large atomic vibration at high temperature [21]. In summary, the rapid crystallization could be explained as follows. Phase-change materials [21,22] present high structural symmetry with rather large distortions. Atoms are randomly distributed in the crystalline phase due to the large atomic displacements at high temperature [23], the crystalline phase is similar to the liquid phase in the sense of "random" atomic distributions. The amorphous phase reflects as a consequence the very isotropic atomic configuration in the liquid state. Finally, the resemblance between the amorphous and the crystalline structures in the atomic distributions at high temperature reduce the barrier between the low states and facilitate the reversible transition through atomic redistribution's and rearrangements. The previously stated macroscopic explanation of the PCM phase crystallization needs to be substantiated by a full understanding of the amorphous structure and the crystallization kinetics. Recently many studies unveiled the details of the transitional mechanisms related to a given PCM [24][25][26][27][28][29] and the nucleus-driven crystallization of amorphous Ge 2 Sb 2 Te 5 [30] and even the pressure induced amorphization of Ge 2 Sb 2 Te 5 [31]. The first clue into the understanding of the crystallization speed was presented by Kolobov & Co. [24] and termed as the "umbrella flip". Using the EXAFS technique, this team was able to reproduce the amorphous state by assuming a tetrahedral atomic arrangement of the Ge atoms, whereas the same atoms have an octahedral-like atomic arrangement in the crystalline phase. (see Fig. [5.3]). This explanation has been refuted by many studies [32]. 4 The partial pair correlation functions: g GaGa (r), g GaSb (r), g GaTe (r), g SbSb (r), g SbTe (r) and g TeTe (r) respectively. Amorphous state (blue lines), liquid state (red line).

The first peak position (FPP), the first minimum position (FMP), and the partial coordina-details) implemented in the package CP2k (model C). It is important to recall that the pseudopotentials employed within Model B and Model C are not the same, being the GTH ones for Model C and the MT for Model B. Figure [6.5] shows the temperature evolution during 300 ps of Models A, B and C. Zero pressure conditions were imposed with careful density calibration and data were collected over a ∼ 25ps of annealing at 300K.

Pair correlation functions and bond angle distributions

Figure [6.6] shows the calculated pair correlation functions of the melt-quenched models, the corresponding FPP and FMP as well the coordination numbers are given in table [6.2].

Overall the shapes of the PCFs for models B and C are different from the one exhibited by model A. The first peak position in g GaGa (r) shifts to higher r values for model B. However, model C gives a better agreement with the values quoted in literature [11,12]. The high intensity of this peak leads to a higher partial coordination number n C GaGa = 0.3703 and a higher number of Ga homopolar bonds. With respect to model A, g B GaGa (r) and g C GaGa (r) show clear second minimums. In the long range region, r > 6 Å, g B GaGa (r) and g C GaGa (r) are more structured than g A GaGa (r). In these new calculations, a given configuration obtained with, say, CPMD is relaxed by using the GPW/SGCPMD approach implemented in CP2k, and viceversa. To this end, three new models are defined as following:

Nearest-neighbor and ring analysis

• Model D: model A (117 atom, melt-quenched in 100 ps with CPMD, detailed above) annealed at 300K for 16ps with CP2k.

• Model E: model B (117 atom, melt-quenched in 300 ps with CPMD, detailed above) annealed at 300K for 16ps with CP2k.

• Model F: model C (299 atom, melt-quenched in 300 ps with CP2k, detailed above) annealed at 300K for 10ps with CPMD.

The partial pair distribution functions thereby obtained are displayed in figure [6.13].