Gheorghe Sorin

Juliana Srour

Fouad El

Haj Hassan

Andrei Postnikov

First and foremost, I would like to thanks my advisors for their guidance and help throughout my thesis -especially prof. Andrei Postnikov whose assistance, patience and support lead to the valuable research included in this dissertation. In fact, he deserves deep and profound gratitudes for many things: notably, for creating a suitable research enviroment, whereby he encouraged me to attend workshops, schools and conferences that enriched my personality in the spirit of doing research. I came to Lorraine University intending to join his group for two primary reasons: the areas of research in which he was involved, and what I had heared about his personality. I have not been disappointed in either regard. The direction of my graduate work has been strongly influenced by the director of our group and our laboratory, Prof. Olivier Pagès, who permanently was ready to help not only on administrative issues related to this thesis, but also by discussing science or life on many occasions. My gratitude goes for his inestimable advises

Résumé

Ce travail expose les résultats de simulations numériques, réalisées dans le cadre de la théorie de la fonctionnelle de la densité (DFT), pour les structures cristallines d'équilibre et pour les propriétés électroniques de quelques semi-conducteurs binaires ou (pseudo)ternaires à base de Cu, In, Ga et Se. Les systèmes étudiés possèdent la même structure à courte portée (environnement tétraédrique des cations et anions) mais diffèrent à longue portée. Les composés binaires (Ga/In)Se, (Ga/In) 2 Se 3 constituent des références importantes dans les diagrammes de phases des systèmes à base de (Cu, In, Se) et (Cu, Ga, Se), au sein desquels figurent les phases potentiellement utiles dans le domaine du photovoltaïque. Le travail comprend deux chapitres d'introduction et trois chapitres exposant des résultats nouveaux.

Le chapitre 1 explique le niveau d'approximation retenu pour la mise en oeuvre de la théorie. Un bref aperçu de la DFT est suivi d'une introduction des méthodes de calcul utilisées (WIEN2k et VASP) et permet d'identifier les paramètres importants du calcul. Une attention particulière est accordée à certains aspects techniques importants pour ce travail : la version PBEsol d'approximation du gradient généralisé (GGA) dans le cadre de la DFT, la variante de Becke-Johnson modifiée (mBJ) pour le potentiel d'échange-corrélation (XC), l'inclusion phénoménologique des interactions dispersives.

Le chapitre 2 expose les problématiques abordées, en liaison avec les structures cristallines des systèmes étudiés, à la lumière des résultats antérieurs, notamment ceux issus de la théorie. La diversité des structures est considérée dans le cadre d'une approche unitaire, couvrant aussi bien Cu(In,Ga)Se 2 , de structure chalcopyrite, que les systèmes lamellaires InSe et GaSe, et enfin les systèmes de lacunes ordonnées, In 2 Se 3 et Ga 2 Se 3 . La construction des supermailles pour les différentes phases ainsi que les relations entre leurs zones de Brillouin sont expliquées.

Le chapitre 3 porte sur les composés Cu(In/Ga)Se 2 de structure chalcopyrite, déjà étudiés par le passé dans la littérature. Dans ce travail, l'attention se porte sur les systèmes Cu(In,Ga)Se 2 mixtes de concentrations In:Ga 1:3, 1:1 et 3:1, étudiés sous forme de structures ordonnées en In et Ga. Les paramètres de maille et longueurs de liaison ont été optimisés pour chaque concentration. D'après les structures de bandes calculées dans l'approximation mBJ, le gap optique augmente progressivement de 0,46 eV dans CuInSe 2 à 1,08 eV dans CuGaSe 2 . Cela reste en deçà des valeurs expérimentales Résumé correspondantes, de l'ordre de 0,8 eV, mais cependant constitue une amélioration par rapport aux estimations basées sur la GGA.

Le chapitre 4 porte sur les semi-conducteurs binaires InSe et GaSe de structure cristalline hexagonale, dont la caractéristique essentielle est la couche double anioncation-cation-anion. Les variétés d'empilement de cette bicouche donnent lieu à des polytypes dont les énergies sont très proches. Les courbes énergie/volume calculées avec les méthodes WIEN2k et VASP, à partir de la même fonctionnelle XC (PBEsol), sont presque identiques. A particulier, nous nous sommes attachés à distinguer les polytypes en fonction de leur stabilité relative. Les différences d'énergies totales entre les polytypes sont très faibles, de l'ordre de 0,5 meV par unité de formule chimique. Les valeurs de bande optique obtenues par l'approximation mBJ ne sont que légèrement sous-estimées par rapport aux valeurs expérimentales.

Le chapitre 5 discute les résultats de calcul pour Ga 2 Se 3 et In 2 Se 3 . Pour chaque composé, trois variantes de structure hexagonale ont été testées : l'une contenant des lacunes cationiques ordonnées en forme de vis et dans les deux autres une couche cationique fait défaut. L'effet de l'ordre particulier des lacunes sur la relaxation atomique et sur la structure de bandes électroniques a été discuté.

List of Tables Introduction

Numerical (computer) simulations of materials have long and glorious history. The simulations done at the microscopic level, i.e. those that "care" about individual atoms and electrons rather than about continuous or effective medium, have proven successful in explaining electron-related properties like magnetism, cyclotron orbits, or X-ray spectra. The calculations done on the basis of very general considerations (i.e., just quantum-mechanical equations), rather than being parametrized for this or other group of materials, are known as "first-principles" calculations. In the last decades, first-principles calculations have seen big success in, for one thing, optimisation of the known crystal structures and even in the prediction of new ones; for the other thing, in the analysis of electron (e.g., optical) excitations. This allowed a broad opening of first-principles calculation techniques towards the problems which arise in the search of new materials and characterizing them under the angle of possible applications for electronics or optics.

The present project started from a general interest in semiconductor-based photovoltaics, and from an ambition, fed by preliminary discussions with my director's colleagues, notably Michael Yakushev (Strathclyde University, Glasgow, Scotland) and Tatiana Kuznetsova (Institute of Metal Physics, Ekaterinburg, Russia), to guess / simulate / characterize, on the basis of first-principles calculations, the phases within the Cu-(In,Ga)-Se system which come about, mostly unintentionally, in the process of synthesis of good materials for photovoltaics. "Good" materials were, and still are, the prototypes like CuInSe 2 and CuGaSe 2 of the chalcopyrite structure, whose "neighboring" phases (in the sense of practical sample preparation), sometimes intentionally added and sometimes annoying yet unavoidable, are copper and indium/gallium selenides. In the course of work, the study of properties of these selenides absorbed more and more effort and eventually led to a reshaping the original plan of work. Namely, the "simple" binary selenides revealed the variety of phases and difficulties in characterizing them by means of calculation, that was recognized as an independent interesting challenge. The focus of the work then shifted onto exploring how the modification of structure brought with itself the modification of the electronic properties. Whereas the shortrange order motives, the tetrahedral environment of cations as well as anions, remained the same in all structures concerned, the ordering on a longer scale was different. For Introduction chalcopyrite-type materials, it resulted in the crystal lattice resembling the zincblende one, only with a variety of cation species distributed over sites. The 1:1 binaries (InSe, GaSe) arranged themselves into layered structures with cation-cation covalent bonds inside and non-saturated (having only three neighbours) outer Se surface. It seemed fascinating that among the many polytypes described in the literature, one particular structure seem to definitely "win" for InSe and a different one -for GaSe. The 2:3 binaries (In 2 Se 3 , Ga 2 Se 3 ) offered a further variety of structures, since they might presumably incarnate either a regular wurtzite, albeit with ordered cation vanancies, or a layered structure with thicker layers. These since long suggested structure models remained rather hypothetical, so their verification by means of first-principles calculations seemed necessary.

Even as the corresponding crystal structures are not very complex, the generation and assessment of meaningful results required accuracy and patience. I wanted to present a study that would result in reliable conclusion; therefore the calculation method initially used was one that is considered as state-of-art of DFT calculations, namely, an all-electron full-potential augmented plane waves method, realized in the WIEN2k code. The requirements for the calculation method, in view of the initial problematics, was the following: (i) to provide an ultimately accurate description of equilibrium structures, at least to overcome the usual underestimation of lattice parameters within the LDA and their frequent overestimation within the GGA; (ii) to get the band structures at least so good as to correctly identify the band gap whenever the material is a semiconductor (and, in the ideal case, to get the gap value about correctly). The tools "routinely" available since certain time to face these challenges are, respectively, the PBEsol flavour of the exchange-correlation (a prescription for the XC functional, that yields the energies and forces and everything that is needed for structure relaxations), and the modified Bethe-Johnson prescription for the XC potential, that provides just the band energies but, consequently, the band gaps. Both are implemented within the WIEN2k. However, the preliminary work done on GaSe and InSe compounds, along with the analysis of information available from previous calculations, revealed that the bonding across the so-called van der Waals gap may come out underestimated in "conventional" DFT calculations, that may result is slight overestimation of the c lattice parameter. It seemed important to include the dispersion interactions into the calculation, that is nowadays possible in a number of computer codes, usually not really from first principles, but on a phenomenological basis. A convenient solution was to use the VASP code which provides this particular option, and has otherwise a reputation of being quite accurate, even if, a priori, not such accurate as the all-electron WIEN2k. As it happened, the larger part of results, namely all those concerning binary systems, was obtained using the both methods, WIEN2k and VASP, in comparison, that helped to reinforce the "consistent" message and get an idea of the "noise". The invaluable assistance of Dr. Michael Badawi in my getting acquaintance with the VASP code, including Introduction numerous discussions concerning details of calculations, was essential for bringing this part of the work to conclusion.

As the work is now done, one can regret that the spin-orbit interaction was completely omitted from the study. Based on the experience of early calculations for some of related systems, one could recognize its relative unimportance in the band structures. However, in view of the differences of total energies between the competing polytypes, which turned to be really tiny, the effect of the spin-orbit on the hierarchy of total energies among the polytypes might have turn out not negligible. The second omission concerns the effect of lattice vibrations on comparing the stability of polytypes. Even as relatively demanding, such calculations could have been in principle done with the VASP code, but was left out because of technical problems and lack of time.

The work consists of two introductory chapters, followed by three chapters that contain the original results.

Chapter 1 defines the place of the present work in the context of previous studies, explaining notably the level of theory applied and its main approximations. A brief overview of the density functional theory (DFT) is followed by a concise introduction into the calculation methods used (WIEN2k and VASP) and is concluded by explaining the essential technical parameters of the calculations and the criteria / routines of their necessary testing. An attention is given to the calculation "features" important for the present work: the PBEsol version of the (meta-)generalized gradient aproximation (GGA), the modified Becke-Johnson (mBJ) approximation for the exchange-correlation (XC) potential within the DFT, the phenomenological inclusion of dispersion interactions.

Chapter 2 explains the problematics concerning the systems under study, their crystal structures and the results of previous works, notably those in theory. The diversity of structures is set into a general frame, passing from Cu(In,Ga)Se 2 with chalcopyrite-type structure to layered InSe and GaSe to ordered-vacancies phases In 2 Se 3 and Ga 2 Se 3 . The construction of unit cells for different phases and the relations between the Brillouin zones are explained, in order to keep the following chapters clear from such technical details.

Chapter 3 deals with CuInSe 2 and CuGaSe 2 compounds of the chalcopyrite structure, which have been already subject to many studies because of their importance for photovoltaic applications. In the present work, the mixed Cu(In,Ga)Se 2 systems at 1:3, 1:1 and 3:1 concentrations of In:Ga have been studied, represented by model supercells that presume an ordered distribution of In and Ga atoms. This allowed to probe the variation of lattice parameters and of the bond length with concentration. The band structures calculated with the mBJ approximation for the XC potential yield the optical gap that underestimates the experimental values by ∼0.8 eV but presents a big "improvement" (in this sense) over "conventional" versions of GGA, remaining at the same time computationally affordable.

Introduction

Chapter 4 addresses the binary InSe and GaSe semiconductors with hexagonal crystal structure, the essential element of which is the anion-cation-cation-anion double layer. An issue of special attention of the work was the discrimination of polytypes according to their expected stability. The energy/volume curves for these systems have been explored on the base of calculations done with WIEN2k and VASP methods. The utmost effort was applied to acquire such differences reliably, by enforcing the necessary accuracy of calculation methods by increasing the corresponding cutoffs. The use of mBJ XC potential results in opening of band gaps, with only slight underestimation of experimental values.

Chapter 5 outlines the calculation results for Ga 2 Se 3 and In 2 Se 3 , comparing, for each of these compounds, three earlier proposed structure models (within the hexagonal setting): that of cation vacancies ordered in screw form and two models with missing cation layer. The impact of the type of vacancies ordering onto the type of relaxation and the electronic band structure are discussed.

The work is terminated by the general conclusion.

Chapter 1

First-principles calculations: overview of the research method

Explaining the level of theory

This work is dedicated to analysis of electronic structure, stability and phase equilibrium of several binary, ternary and mixed semiconductor compounds. The importance of these systems and the related problematics will be discussed in Chapter 2. Here, we briefly discuss the methods by which the systems are studied. The electronic structure implies quantum-mechanical approach, the calculation of energy bands and thence derived total energies. At earlier stages of studies on semiconductors, band structures were derived from model Hamiltonians, the parameters of which were extracted from some crucial comparisons with experiments. Gradually, electronic structure calculations referred to as "first-principles", or ab initio ones, became feasible and gained in accuracy. "First principles" means that these calculations do not depend on parameters tuned specifically for a particular system, but are performed on the basis of quite general foundations of quantum mechanics. In the ideal case, the ab initio perspective implies that it suffices to pick the constituent atoms in a given relation, and then their interactions, resulting energy levels, potential energies, forces acting on atoms, and, consequently, the equilibrium structures can be predicted fully within the framework of theory, driven by the equations of the quantum mechanics.

On the way of practical realisation, there are different possible levels of accuracy that need to be explained here, and the level used in the work specified.

Speaking of quantum mechanics, we recall the Schrödinger equation and its wave function, dependent on the coordinates of all intearcting particles, that are, from the point of view of condensed matter problematics, all electrons and nuclei of the system. The analytical solution is probably limited by the case of (one-electron) H + 2 molecule [see, e.g., [START_REF] Grivet | The hydrogen molecular ion revisited[END_REF] whereas already the (two-electron) helium atom can only be Chapter 1 treated numerically [START_REF] Tanner | The theory of two-electron atoms: between ground state and complete fragmentation[END_REF]. The first level of broadly acceptable approximation is that of [START_REF] Born | Zur Quantentheorie der Molekeln[END_REF], or the adiabatic one, that allows to separate the movements of electrons and nuclei on the grounds of big difference in their masses. 1 For all practical purposes, the use of the Born -Oppenheimer approximation amounts to the solution of Schrödinger equation for electrons in the electrostatic field of fixed nuclei. Further on, the forces induced on nuclei by other nuclei and by electrons can be evaluated and used to correct the nuclei positions, eventually leading to structure optimisation.

Speaking henceforth of electronic problem only, one faces the difficulty that a search for the many-particle wave function, even by numerical means and already for a very limited number of electrons, becomes a tremendously complicated problem. Historically, a milestone was set by the Hartree -Fock (HF) method, or formalism [see, e.g., Echenique and Alonso, 2007, for a review], that suggested an ansatz for the many-particle wave function constructed as a Slater determinant of one-electron wave functions. By itself, the HF is a mean-field method and hence not particularly accurate one, albeit still largely in use for practical calculations. However, the HF method is important because it sets a stage for gradual improvement of aproximations for many-electron wave function, taken as combinations of Slater determinants. Such developments make a basis for different advances in quantum chemistry; they reached big sophistication and became very accurate in relation to molecules. On the contrary, the generalisation of such methods to extended systems, e.g., crystals, faces serious problems.

A different way to numerically construct a many-electron wave function is offered by Quantum Monte Carlo (QMC) methods -see, e.g., [START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF] for a review.

It turns out that the knowledge of the wave function is not necessary for a number of practical applications. Notably, the total energy, the central property in performing structure optimisations and thermodynamic analysis, can (in principle) be unambiguously expressed, as it is a functional, of electron density, itself an important observable property. This is true if applies to ground state of system, hence a priori excludes excitations, but remains useful for a great number of practical applications, including search for equilibrium structures, study of phase diagrams, lattice dynamics. The link total energy -density is the core of the density functional theory (DFT) as it took shape since the pioneering works of Walter Kohn, Pierre Hohenberg and Lu Sham [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]].

Elements of the Density Functional Theory

1.2 Elements of the Density Functional Theory

General context

The DFT being the essential theory base of the calculation methods employed in the present work, it deserves a more detailed introduction although a vast literature pool exist in the subject, offering different levels of immersion [see, e.g., [START_REF] Burke | The ABC of DFT[END_REF][START_REF] Capelle | A bird's-eye view of density-functional theory[END_REF][START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]. The historical roots of the DFT are the Thomas -Fermi theory [see [START_REF] Spruch | Pedagogic notes on Thomas-Fermi theory (and on some improvements): atoms, stars, and the stability of bulk matter[END_REF], for a review] where the charge density for the first time enters as an esential parameter in the quantum-mechanical description of the electron gas, and the HF theory, starting from which the extensions have been undertaken, by means of perturbation theory, towards microscopic description of manyparticle systems.

In a nutshell, the work by [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] suggests a practical scheme how the total energy and the density can be calculated. Namely, the density sought for, ρ(r), is cast as that due to non-interacting quasiparticles, and hence constructed from their single-particle wave functions,

ρ(r) = N i (occ.) |φ i (r)| 2 , (1.1)
where the summation is done over the lowest N (= number of electrons) states, namely the occupied ones according to the aufbau principle, in the spectrum of the system. The (fictituous) functions φ i (r) and their corresponding eigenvalues ǫ i , needed to sort out the solutions and to apply the aufbau principle in Eq. (1.1), are given by the Kohn -Sham (KS) equations (in atomic units, = electron mass = elementary charge = 1):

- 1 2 ∇ 2 - nuclei µ Z µ |R µ -r| + ρ(r ′ ) |r -r ′ | dr ′ + v xc (r) φ i (r) = ǫ i φ i (r) . (1.2)
These equations are very similar to the HF ones, only that the (nonlocal) exchange term in HF is replaced by (local) exchange-correlation (XC) potential v xc (r) acting at each orbital φ i (r). Being local (i.e., the same for any orbital at which it acts), the XC potential makes the solution of KS equations technically less demanding than that of non-local HF equations. However, the v xc (r) is in fact a functional of density and may thus depend on r in a quite complicated way. Even worse, the shape of v xc (r) is not exactly known, even if a number of its special and asymptotic cases have been studied. Essentially, v xc (r) incorporates all effects due to "multiparticularity" which enter on top of kinetic energy, the first term on the left of Eq. (1.2), and the mean-field electrostatic interaction with the whole electron density, the third term in the left. v xc (r), a priori Chapter 1 much smaller than each of two mentioned terms, corrects for self-interaction (the absence of the electron concerned in the charge density with which it interacts) and takes into account how do the electrons participating in the charge density avoid each other, due to the Pauli principle (that is, exchange, already present in the HF formalism), and due to particle-particle (beyond the mean-field) electrostatic repulsion (that is correlation, absent in the HF). In fact, all the complexity of interactions described by v xc (r) can be understood as a pure electrostatic interaction between the charge density and the exchange-correlation hole, the residual density expelled from the vicinity of each individual electron. Some important conclusions concerning the spatial extension and the shape of the XC hole can be done on the basis of electron gas models, whereas more specific calculations for reference cases beyond simple limiting cases can be obtained from numerical QMC calculations. Such, [START_REF] Ceperley | Ground state of the electron gas by a stochastic method[END_REF] scanned the phase diagram of electron liquid throughout a broad range of densities and phases (ferromagnetic of antiferromagnetic gas, Wigner crystal, etc.) This simulation laid foundation for a number of practically useful fitting of the XC potential. In fact, a given model describing the XC hole depending on electron density can be straightforwardly transferred onto the resulting v xc (r).

Practical recipes for the XC potential: LDA and GGA

The XC potential often appears in discussions as something mysterious, that is known to exist but whose exact shape escapes a precise specification. Certain cornerstones can be indicated to facilitate the orientation on this terrain. First, it is convenient to perform a separation between the exchange v x and the correlation v c parts, the first being fully grasped by the HF approximation, and the second remaining fully beyond it. The exchange part can be, in principle, exactly calculated for a system of electrons whose one-particle functions φ i (r) are known. The exact result for the homogeneous electron gas of a uniform density ρ, following from the integration of the HF exchange over all occupied free-electron states, is [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]. This result is a starting point of generalisation over the case of slowly varying density, where the constant prefactor was subject to ad hoc adjustment, as in the X α method [START_REF] Schwarz | Approximate numerical Hartree-Fock method for molecular calculations[END_REF], and eventually gave place to more "universal" prescriptions containing some smooth function of ρ(r) in the prefactor, that would allow a reasonable interpolation between different reference cases (e.g., those of fully spin-polarized and spin-compensated electron gas). The correlation part is more difficult to extract; the corresponding benchmark results for construction of working formulas were initially gained from many-body perturbation theory expansions,

v x = - 1 π [3π 2 ρ] 1/3 . (1.3) -see Eq. (2.19) of

Elements of the Density Functional Theory

and then, starting from the seminal work by [START_REF] Ceperley | Ground state of the electron gas by a stochastic method[END_REF], -from QMC simulations. The generic name local density approximation (LDA) covers all realisations in which the quite general form for the XC energy as some functional of density

E xc [ρ] = dr ǫ xc (r) ρ(r) ; v xc (r) = δE xc δρ(r) (1.4)
(without loss of generality, the XC energy density ǫ xc , itself a functional of density, is weighted over space according to the distribution of ρ) is replaced by the simplified form

E LDA xc = dr ǫ xc ρ(r) ρ(r) , (1.5)
in which the XC energy density ǫ xc (r) is presumed to be just a function of ρ in the r point only. More precisely, an explicit dependence on two spin-resolved components of charge density, ρ ↑ (r) and ρ ↓ (r), may be allowed, that is referred to as local spin density approximation (LSDA). The practical parametrisations are often expressed in terms of the Seitz radius r S = [3/(4πρ)] 1/3 i.e., the radius of sphere containing one electron at the given density ρ = ρ ↑ +ρ ↓ (a standard measure of "electron dilution" in DFT works), and the relative spin polarisation ζ = (ρ ↑ -ρ ↓ )/ρ, that varies, in each given r, between 0 and 1:

E LSD xc [ρ ↑ , ρ ↓ ] = dr ǫ xc (r S , ζ) ρ(r) . (1.6)
Historically important examples of the LDA parametrizations are those of Hedin-Lundqvist [START_REF] Hedin | Explicit local exchange-correlation potentials[END_REF]], Perdew-Wang 92 [START_REF] Perdew | Accurate and simple analytic representation of the electron-gas correlation energy[END_REF]] and a number of others in between. In practical terms, calculation packages usually offer certain variety of these options; however, the results of calculations with different LDA parametrisations come out quite similar. Under "results", we mean here the Kohn-Sham band structures ǫ(k) calculated for a given arrangement of atoms, or geometric structure (unit cell shape, bond lengths) obtained from the minimisation of the total energy if the calculation code allows this. The reasons of unexpectedly good performance of LDA, in spite of its quite drastic nature, are discussed in many works on DFT, e.g., [START_REF] Jones | The density functional formalism, its applications and prospects[END_REF]. A notorious drawback of LDA is its known "overbinding", that is, a too short (typically by just few per cent) prediction of equilibrium lattice parameters and interatomic distances, and a bit too large estimations of bonding energies. An improvement can be achieved by taking into account not just density at the point r at which ǫ xc (r) must be defined, but also the gradient of density, i.e., ǫ xc (r) is casted as ǫ xc ρ(r),∇ρ(r) . Such an extension is usually referred to as nonlocal corrections, to make a distinction from a genuine non-local functional form ǫ xc [ρ(r)] that would depend on the shape of charge density over some extended spatial region. The practical realisation is more sophisticated than just formally adding a term of the Taylor expansion. In fact Chapter 1 the additional parametrisation permits to satisfy at least some features of the correct formal behavior expected from the XC potential in certain limiting cases, that was not granted within the LDA. For instance, the correct asymptotic of v xc (r) on removal an electron from the rest of system must be Coulombic, v xc ∼ -1/r for r → ∞, dominated by the attraction towards the net charge of the system left behind, whereas the LDA XC-potential falls down as ∼ -[ρ(r)] 1/3 and hence exponentially, i.e., an electron is too weakly bound. This wrong asymptotic is the origin of some deficiencies of LDA; an enlightening discussion was offered, e.g., by van Leeuwen and Baerends [1994]. A variety of forms that specify the functional dependence of the XC energy density in terms of exactly four arguments, i.e., the spin-split densities and gradients, subject to certain "physical" constraints, fall into the category of generalized gradient approximation (GGA), discussed at length by Perdew et al. [1992]:

E GGA xc = drf (ρ ↑ , ρ ↓ ,∇ρ ↑ , ∇ρ ↓ ) , (1.7)
Here again as with the LDA, a number of convenient parametrisations has been proposed and tested. We single out, for the following reference, the Perdew-Burke-Ernzerhof (PBE) scheme [Perdew et al., 1996] which became a de facto standard of GGA calculations, being cited in over 42400 articles. The explicit form of the XC energy according to the PBE-GGA is

E PBE xc [ρ ↑ , ρ ↓ ] = drρ(r) ǫ unif x (ρ(r)) F xc (r S , ζ, s) , (1.8)
where F xc (r S , ζ, s) is the enhancement factor of (non-local) exchange-correlation over the (local) exchange density of the uniform density distribution (that for which the LDA is exactly valid):

ǫ unif x = - 3 4 3 π 1/3 ρ 4/3 , (1.9)
the functional derivative of which over ρ returns (1.3). The enhancement factor F depends, along with earlier defined Seitz radius r S and spin polarisation ζ, also on the "non-locality" s that a dimensionless measure of the density gradient at given point r:

s = |∇ρ| 2k F ρ , with the local Fermi radius k F = (3π 2 ρ) 1/3 = 1.91916/r S .
(1.10)

The "non-locality" may vary from s = 0 (hence passing into the LDA limit) to s 3 "for real systems", citing Perdew et al. [1996]. The parametrisation of F xc (r S , ζ, s) is discussed by Perdew et al. [1996], and the shape of this function in dependence of all three its arguments is depicted in Fig. 1 [1996] that the small-s limit of the enhancement function for exchange must be F x (s → 0) = 1 + µs 2 for all GGA schemes, so that technical discussions concentrate on the value of the coefficient µ and on the possible choices to conveniently specify F x (s) throughout the broad s range, e.g., Perdew et al. [1996] suggested

F x (s) = 1 + κ - κ 1 + µs 2 /κ , (1.11)
with the constants µ = 0.21951 and κ = 0.804. In fact some of parameters are defined by inequalities rather than by exact fit, that leaves a certain freedom in parametrization, that was used to enforce possibly best agreement of the calculation results for a large number of chosen benchmark systems.

A modified parametrisation of F xc (s) has been proposed by Perdew et al. [2008b], who aimed at obtaining more accurate results (equilibrium latice constants; atomization energies) for bulk solids and surfaces, raher than for isolated atoms to which the PBE case was tuned. Their formula became largely used under the name PBEsol. As is stated in this work, "PBEsol should improve most surface energies over LSDA, whereas PBE worsens them", and further on "The nonlocality or s dependence of GGA exchange is diminished from PBE to PBEsol, making the latter somewhat closer to LSDA". On the web site of the Kieron Burke's group, one can find documented Fortran codes of the subroutines to calculate the XC potentials and energies for the PBE (http:// dft.uci.edu/pubs/PBE.asc) and PBEsol (http://dft.uci.edu/pubs/PBEsol.html) algorithms. The following fragment refers to the definitions covered by the formulae (1.9), (1.8) and (1.11) above. The difference between PBE and PBEsol is only in the parameter values applied in one or the other case. ---------------------------------------------------------------------...

... c---------------------------------------------------------------------- c Formulas: c e_x[unif]=ax*rho^(4/3) [LDA] c ax = -0.75*(3/pi)^(1/3) c e_x[PBE]=e_x[unif]*FxPBE(s) c FxPBE(s)=1+uk-uk/(1+ul*s*s) [a](13) c uk, ul defined after [a](13) c-
An excellent discussion on the backgrounds and implementation details of the GGA can be found in the Chapter 18 of Burke and friends [2007].

Kohn -Sham eigenvalues and band gaps

The results of DFT calculation are not just total energy or charge density, but also the KS eigenvalues. If the system under study is dielectric or semiconductor, its important characteristic is the band gap (and also the characterisation of the band gap as either direct or indirect one). The latter is related to energy dispersion, see section 1.3.3. In Chapter 1 experiment, the band gap can be extracted from optical measurements, and the band structure recovered from angle-resolved spectroscopy. In the community that deals with DFT calculations, there is a tendency to compare the "experimental" (spectroscopic) band structures with the dispersion of KS eigenvalues, and to compare the optical band gap with the energy difference between the highest occupied and the lowest unoccupied KS orbitals. This tendency can be explained by a certain similarity between the KS and HF context and the shape of resulting equations. In the HF formalism, the Koopmans' theorem holds that relates the energy of a given HF orbital to the difference of total energies of the system in two cases, with the orbital in question being occupied or vacant. In the KS formalism, the eigenvalues do not have this meaning. In principle, one should take into account that the DFT was largely derived keeping in mind the ground-state properties, whereas the optical band gap is something characterising the electron excitation. However, the relation between the two exist and was elucidated in a number of works, e.g., quite nicely in a review by [START_REF] Capelle | A bird's-eye view of density-functional theory[END_REF]. The physical meaning of KS eigenvalues is given by the Janak's theorem [START_REF] Janak | Proof that ∂E/∂n i = ǫ i in density-functional theory[END_REF] as the derivative of the total DFT energy, E, with respect to the occupation of a given orbital n i :

ǫ i = ∂E ∂n i (1.12)
In DFT, the total energy E can be generalized over non-integer number of electrons N . [START_REF] Perdew | Density-functional theory for fractional particle number: Derivative discontinuities of the energy[END_REF] have shown that E(N ) is a piecewise linear function and hence not analytical, having kinks (slope discontinuities) at integer N values. This discontinuity at the factual number of electrons N enters the theoretical expression for the optical gap ∆,

∆ = δE[ρ] δρ(r) N +δ - δE[ρ] δρ(r) N -δ (1.13)
that in its turn can be attributed to two contributions to the total energy from the kinetic energy T [ρ] and the XC energy:

∆ = δT [ρ] δρ(r) N +δ - δT [ρ] δρ(r) N -δ ∆ KS = ǫ N +1 -ǫ N + δE xc [ρ] δρ(r) N +δ - δE xc [ρ] δρ(r) N -δ ∆xc = v + xc (r) -v - xc (r)
.

(1.14)

The "Kohn-Sham gap" ∆ KS immediately follows from a solution of the KS equations; it is "underestimated" in the sense that the second term, ∆ xc , that should be provided by the "exact" XC potential, is not a priori included in the practical realisations used (as those discussed in the previous section). Any XC potential designed to "correct" the band gap on the physical basis and not just empirically should somehow imitate the above mentioned slope discontinuity of E xc (N ).

Elements of the Density Functional Theory

The version of v xc implemented to this end in the WIEN2k code (used in this work to calculate band structures and to discuss the band gaps) was suggested by [START_REF] Tran | Accurate band gaps of semiconductors and insulators with a semilocal exchange-correlation potential[END_REF] under the name "modified Becke -Johnson" (mBJ), because this prescription is based on the paper by [START_REF] Becke | A simple effective potential for exchange[END_REF]. The overview of the performance and the accuracy of the mBJ potential based on early examples of its use was given by [START_REF] Koller | Merits and limits of the modified Becke-Johnson exchange potential[END_REF].

The kinked stepwise behaviour of E(N ) indirectly manifests itself in the shell structure of atoms; Figs. 1,2 of [START_REF] Becke | A simple effective potential for exchange[END_REF] show the "waves" in the radial profile of "accurate" exchange potentials which get very much damped in the LDA. It is technically difficult to re-construct such oscillations having at one's disposal just the charge density, which is quite smooth in space yet varies over orders of magnitude (although, in the spirit of the DFT basics, the change density must somehow contain all the necessary information). The idea of Becke and Johnson was to use the kinetic energy density as a robust and efficient tool to inspect where the undulations of the exchange potential need to be added in order to reconstruct the desired atomic-shell related structure. The working formula must be quite general and not related to particular atoms of electron shells. [START_REF] Becke | A simple effective potential for exchange[END_REF] suggest to use the (spin σ -resolved) kinetic energy density τ σ ,

τ σ (r) = 1 2 Nσ i(occ.) |∇φ iσ (r)| 2 , (1.15)
normalise to the density itself, and add it, with a prefactor chosen to satisfy the uniform electron gas limit, to the "conventional" reference exchange potential v ref.

x (that of Slater, in their paper):

v BJ xσ (r) = v ref. xσ (r) + 1 π 5 12 2τ σ (r) ρ σ (r) . (1.16)
The modification of this scheme due to Tran and Blaha [2009] was that they suggested to use, as a reference exchange potential, that of [START_REF] Becke | Exchange holes in inhomogeneous systems: A coordinate-space model[END_REF], expressible via a solution of a non-linear equation depending on the kinetic energy density, charge (spin) density, the density gradient and Laplacian. Moreover they introduced a variable prefactor (depending on the average value of |∇ρ|/ρ) that affects the relative weights of two terms in the right-hand side of Eq. (1.16). The remaining freedom in fixing this prefactor was used to optimise, on the average, the calculated values of optical gaps over large number of insulators and semiconductors as compared to corresponding experimental values. Like in the case of practical schemes for LDA or GGA, such "fixing" is done once for all systems and "works" reasonably for very different systems with different band gaps.

Chapter 1

The generic name for XC potentials which make use, along with (local) density and its gradient, also the kinetic energy density, is meta-GGA [START_REF] Tao | Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids[END_REF]. The mBJ scheme addresses only the exchange potential; the effect of correlations can be added within the LDA but, according to [START_REF] Tran | Accurate band gaps of semiconductors and insulators with a semilocal exchange-correlation potential[END_REF], has only small effect. It should be noted that the use of mBJ functional primarily makes sense to get the "improved" band structures, in view of their comparison with experiment. As F. Tran put it in one of his comments to the users' mail archive related to the WIEN2k code: 1 The mBJ potential should be used only for the calculation of the electronic and magnetic properties, but not for quantities which involve the total energy like the elastic properties. MBJ is only a potential and there is no associated functional (LDA was arbitrarily chosen as the energy funtional for MBJ, option 28).

Semiempirical inclusion of van der Waals interaction

The III-VI semiconductors addressed in the present work have double-layer structure; among the interactions between adjacent double layers that keeps the crystal together, an important contribution comes from the van der Waals (vdW) interactions. More precisely, these are London dispersion forces, "...arising from quantum-induced instantaneous polarization multipoles in molecules. They can therefore act between molecules without permanent multipole moments" (Wikipedia). Indeed, the double layers are electrically neutral, do not possess dipole moments, and there is hardly a covalent bonding across the gap between the layers (the so-called vdW gap), because the bonds are saturated within the double layer. The vdW interactions exist between all atoms, but remain negligible against covalent or ionic bonds if the latter are present. However, the omission of vdW interactions may result in an appreciable error in estimating the equilibrium distances between neutral atoms or molecular fragments.

A "conventional" DFT calculation (of the level of LDA or GGA) a priori does not contain the physical mechanisms that would account for the dispersion interaction. The latter was added in the VASP code (discussed below) semiempirically, as described by [START_REF] Bučko | Improved description of the structure of molecular and layered crystals: Ab initio DFT calculations with van der Waals corrections[END_REF]. The formalism realises the "DFT-D2" approach by [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF], where the dispersion interaction is added to the "conventional DFT" part

E DFT+disp = E DFT -s 6 i,j C ij 6 f (|r i -r j |) |r i -r j | 6 , (1.17)
with its r -6 dependence on the distance between the atoms at r i and r j is explicitly inserted, the distance-dependent analogue of the Fermi function

f (r ij ) = 1 1 + exp[-d( r ij R ij -1)] (1.18)
1 http://www.mail-archive.com/wien%40zeus.theochem.tuwien.ac.at/msg08020.html 1.3. Practical solution of Kohn -Sham equations effectively cuts off the interaction beyond the distance R ij that is taken as a sum of vdW radii for atoms i and j (tabulated); the atom-dependent pre-factor C ij 6 = C i 6 C j 6 also comes together from tabulated values; the general prefactor s 6 and the smoothness parameter in the Fermi function d are taken as suggested by [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF]. We note that the Grimme correction only affects the total energies (and forces, in the process of performing the structure optimisation). There is no band structure or density of states associated with the vdW interaction.

1.3 Practical solution of Kohn -Sham equations 1.3.1 All-electron treatment vs exclusion of core states

Equations (1.
2), taken together with the definition of density Eq. (1.1) and taking into account a specific expression for v ex (r) as a functional of ρ(r), make a coupled integrodifferential system; their "straightforward" numerical solution by "standard" methods developped for differential equations would be quite cumbersome. At this stage, further approximations need to be taken, for the sake of practical feasibility. The first important "bifurcation" affects the number of states to be taken into account in the course of solving the KS system. The alternatives are, either taking all possible electronic states into account, including those stemming from the most bound ones on each constituent atom, from the 1s electrons upwards, or restricting oneself to some higher, "chemically relevant" electronic states only. There are arguments in favour of each choice, and the decision taken would largely depend on the context and the objective of study. The electron band structure, optical transitions and spectroscopy, magnetic moments, chemical bonding and the distribution of electron density in the interatomic space are all properties which are dominated by valence electrons from higher (valence) shells of constituent atoms, which enter in some interaction. To be specific, for the elements studied in the present work, these states are 4s, 3d and 4p for gallium (13 electrons, distributed over 1×s + 5×d + 3×p → 9 states in an isolated atom), the same amount for its homologue indium, 16 electrons (three 4p electrons more than for Ga) in case of selenium, 11 valence electrons stemming from five 3d and one 4s valence states for copper atom. In comparison, including all electrons would imply 31 per gallium atom, 34 per selenium, 29 per copper and 49 per indium, and the number of electronic states (i.e., the Kohn -Sham orbitals) over which these electrons need to be accommodated, would increase considerably, with the corresponding impact onto the total size of the Kohn -Sham system. The "neglection" of deeper, core electrons presumes that these latter are indeed present in the atoms, localize themselves in the inner 1s, 2s etc. shells, but their effect amounts to just screening the nuclear charge felt by the valence electrons which only do enter into the interaction described by the KS equations. The arguments Chapter 1 in favour of considering the deeper core states along with the outer valence ones is that the shape of the core states, their exact energy position and hence ionisation energy are somehow, even if slightly, affected by the exact distribution of the valence states, so it will affect the analysis of such core states -related experimental properties as Xray spectroscopy, or the hyperfine interactions (of the 1s electrons with the nuclear spin). Moreover, an important in many practical applications, the correct inclusion of core states into the calculation has an effect on the total energy of system, that is often important to make correct predictions concerning competition (relative stability) of different crystallographic or magnetic phases. Therefore the choice between allelectron and "valence-states-only" methods is a matter of compromise between the affordable calculation complexity and the ultimate precision required.

In the present work, two methods are used, in many cases applied to the same system, that allows to draw some conclusions about their comparison. The first method, realized in the WIEN2k package [START_REF] Blaha | WIEN2k, an augmented plane wave + local orbitals program for calculating crystal properties (Karlheinz Schwarz[END_REF], see also http://www.wien2k.at], falls into the category of all-electron ones; we'll come to discuss the details of its realisation in a minute. The second method, realized in the VASP package (http://www.vasp.at), effectively neglects the core states and treats the valence electrons only. Now, the ways to get rid of the core states are different, and not trivial. Usually, it is done within the frame of pseudopotential approach, that re-formulates the potential felt by an electronic state in question, introducing (into a Schrödinger equation, or a Kohn -Sham equation) a term that is orbital-specific and hence non-local. The construction and testing of pseudopotentials is a vast area that is fully beyond the scope of the present work. A database of professionally prepared and tested pseudopotentials is provided along with the VASP code. It should be mentioned that the pseudopotential method does generally provide, for a given atom, a solution (one-electron wavefunction) without radial nodes, that is, yielding a wrong node structure. For example, the radial 4p-like solution for Ga would depart from zero at the origin (nuclei) and fall down beyond passing a single maximum, whereas the genuine 4p function in an atom would have 2 nodes (change sign two times before asymptotically falling down at large distance). This deficiency is, from one hand, physical (since there are no "previous" states of the same orbital quantum number to which the valence state in question needs to be orthogonal); from the other hand, it is intentionally "build in", driven by the desire to avoid valence functions to be strongly fluctuating around the atom cores, that would permit to decompose them into a plane-wave basis with "not-too-high" cutoff (see below). In fact, one of principal objectives in the historical context of developping the pseudopotential methods (from the model pseudopotentials in the textbooks all the way to modern ab initio pseudopotentials) was to get rid of these annoying fluctuations and to advance, as far as possible, within the convenient plane-wave picture. An obvious deficiency of having a wrong node structure is tolerated under understanding 1.3. Practical solution of Kohn -Sham equations that this deficiency is confined to some region relatively close to the nucleus,1 whereas the chemical bonding, band structure etc. are essentially mostly shaped by overlaps of outer "tails" of atom-centered wave functions, which are fairly enough represented by a pseudopotential calculation. Important conclusions to be drawn from these very general considerations concerning the use of pseudopotentials are that (i) the electron energy bands stemming from the states in atoms more deep than those explicitly treated as "valence states" are not present whatsoever, nor do they contribute to the total energies;

(ii) the charge density plots created in a pseudopotential calculation may be reasonable and accurate in the interatomic region, but do not reflect the node structure of genuine wave functions in the vicinity of nuclei, and hence are completely wrong in these areas.

The pseudopotential paradigm is not the only way to get rid of core states mixing up with the valence ones in the calculation. If fierce spatial fluctuations of valence electron states do not pose a problem by themselves (say if these functions do not need to be numerically represented by plane waves), they can be left as they are, with correct nodal structure and resulting in a "correct" charge density (see above). Simultaneously, the deeper states identified as core ones can be simply excluded from the coupled Kohn -Sham equations and solved separately, as for a single electron in a given potential well (as for an isolated atom, or -better -in a resuting (overlapped) potential at the given site in crystal. Such treatment of core states, done "once and for all", is typically referred to as frozen core approach. Yet another -and more modern -technique, effectively combining the advantages of pseudopotential and frozen core ones, is the projector augmented wave (PAW) method elaborated by [START_REF] Blöchl | Projector augmented-wave method[END_REF].

In a nutshell, in place of straightforward pseudofunctions, smooth around the core and "knowing nothing" about their would-be-correct nodal structure, the method employs the functions "flattened" by a special linear transformation, taking into account the shape of deeper orbitals to which these valence ones have to be orthogonal. In this way, a back and forth transformation between "genuine" (fluctuating) and softened / nodeless (projected) functions is defined, that uses the information about the core states as a transformation / projection "key". In this sense, the PAW method is comparable with, but more sophisticated and flexible than, the frozen core technique. The PAW is implemented in VASP and covered by the documentation to the code; moreover a fair Wikipedia entry exists 2 that indicates a number of useful references. No optimisation Chapter 1

or tuning in this sense was undertaken in the present work and default VASP settings were used, therefore this subject won't be specifically covered further on.

To conclude with the issue of separation into the valence and core states, it can be noted that sometimes this discrimination is not easy, if the states in question are localized, but not quite so strongly as the yet deeper lying states. In this case one speaks of semicore states. Among the elements appearing in the present work, the Se3d or Cu3p states can be given as examples. In a solid, such states form a split-off band that does not mix up with other electronic states; yet the related wavefunctions centered at different such sites in crystal do a bit overlap, so that a tiny dispersion (hence bandwidth) of the corresponding energy band can be noted. The inclusion of such states as the valence ones, always possible in all-electron methods, is usually important for obtaining more accurate results for the total energy; this way it was done in the present WIEN2k calculations. In pseudopotential calculations, as indicated above, it is usually possible to include only one state of the given orbital moment (l value) in the valence band, hence a choice must be done whether it should be, e.g. rather Se3d (essentially occupied) or Se4d (essentially empty). This choice can be done ad hoc, driven by the priorities of calculation and the availability / quality (accuracy) of corresponding pseudopotentials. In the VASP calculations done in the present work, the following valence configurations of elements were used (the states, and electron numbers from each atom, that contribute to the solution of the Kohn -Sham equations): Ga3d 10 4s 2 4p 1 , In4d 10 5s 2 5p 1 , Se4s 2 4p 4 .

Basis functions; calculations with WIEN2k and VASP

As was already mentioned at the beginning of the previous subsection, the straightforward solution of the KS system (1.2) by methods applicable for differential equations is hardly practical. The alternatives could be, (i) a discretization of all functions on a regular, or adaptive [see, e.g., [START_REF] Gygi | Electronic-structure calculations in adaptive coordinates[END_REF] spatial mesh throughout the periodic unit cell in case of a crystal, or throughout the "volume" (e.g., a "simulation box") in case of an isolated molecule, or cluster; (ii) the decomposition of KS functions searched for over a system of either fixed, or adjustable, basis functions, and the reduction of the initial integro-differential system to a system of algebraic equations, solvable by means of diagonalization. To specify the approach, we cast the Kohn -Sham system in a compact form incorporating all the terms in the brackets at the left of Eq. (1.2) into the Kohn -Sham Hamiltionian Ĥ, whatever the specific form of the XC term, resulting in Ĥφ i = ǫφ i . Further on, the searched for KS functions are expanded over the basis system |χ α with coefficients C iα . With this, Eq. (1.2) takes the shape

α C iα Ĥ|χ α = ǫ i α C iα |χ α (1.19)

Practical solution of Kohn -Sham equations

that, multiplying each term on the left by

χ β |, yields ∀β, α χ β | Ĥ|χ α ≡ H βα C iα = ǫ i α χ β |χ α ≡ S βα C iα , or HC = ǫ SC , (1.20)
that is, a generalized diagonalisation problem with right-side matrix S, eigenvalues ǫ and eigenvectors C. Knowing the latter permits to recover the KS functions in the original (spatial) representation, taking the basis functions with the now specified coefficients.

With respect to the above decomposition, several remarks can be done.

• There is no basic difference between the discretization on the grid, i.e., the option (i) above, and the decomposition into basis functions, i.e., option (ii), because the basis, as a valid choice, may consist of "primitives" equal to 1 at a given grid point and 0 elsewhere. Then the expansion φ i = α C iα χ α is simply a piecewise (linear, or more sophisticated, depending on "primitives") interpolation of φ i on the grid, and the whole machinery of Eq. (1.20) applies.

• The basis does not need to be orthogonal (that might be anyway difficult to impose, apart from special cases, like plane waves); the presence of the overlap matrix S βα takes this into account.

• The size of the basis ought to be at least as large as the number of occupied states φ i need to be calculated, but in reality (unless the basis is "perfect") it must be much larger, in order to provide the necessary flexibility in adjusting the solutions φ i . As the KS equations are variational, the "improvement" (in practical terms -a systematic increase) of the basis system allows to "improve" (i.e., to obtain lower total energy, otherwise enhance the precision of) the solutions and all related propertiesband energies, charge density, etc.

• The charge density follows from Eq. (1.1); expressed with the help of the basis functions χ α , it takes the form:

ρ(r) = N i (occ.) |φ i | 2 = α,β χ β | N i (occ.) C iβ C iα ≡ ρ βα |χ α , (1.21)
which defines the density matrix ρ αβ in the representation by the basis system |χ . Quite common choices of basis functions for practical electronic structure calculations are:

(1) plane waves (predominantly in combination with pseudopotentials or PAW technique, when the functions φ i do not exhibit strong spatial fluctuations);

(2) atom-centered more or less localized functions that possess either "reasonable" physical shape (Slater-type orbitals, numerical solutions of some reference systems, Chapter 1 etc.), or "convenient" mathematical properties (Gaussian-type orbitals that permit an analytical evaluation of multicenter integrals, e.g., H βα and S βα in Eq. (1.20) ).

In the VASP code, the planewave basis is used. It should be noted that the use of planewaves as basis implies the exact periodicity of the system under study. In the present work, this was naturally given, since all systems are crystalline and characterised by compact enough translation vectors. To apply a planewave code for the treatment of a molecule, it is necessary to place the latter into an artificial "simulation box", imposing a spurious periodicity (and hence "contaminating" the results by spurious interaction of translated replicas of the system under study). The periodicity of the basis functions at the lengths of the simulation cell (or, the primitive cell in a crystal) is automatically guaranteed if the plane waves are numbered by the vectors of reciprocal lattice G, since, for any lattice translation R, holds e iG(R+r) = e iGr . The planewave basis functions can be systematically generated and enhanced by the largest G magnitude taken into account, that is, the wavevector of the fastest fluctuating plane wave. In practice, a cutoff is set as the radius of sphere in the reciprocal space, G max , within which all the vectors G are used for the construction of the basis. Such cutoff is usually expressed in units of reciprocal length, or, using the relation

E cut = G 2 max 2m , (1.22)
-in units of energy. The actual number of the planewaves used scales with the volume of the periodic cell, Ω, as

N PW ∼ ΩE 3/2
cut . An important element of the WIEN2k code is the muffin-tin geometry, a remnant of the historically used muffin-tin approximation (the potential being spherically symmetric within atom-centered non-overlapping spheres and flat in the interstitial region). The muffin-tin (MT) geometry preserves this division of space, although nowadays the intra-sphere potential includes non-spherical terms (it is expanded over spherical harmonics around each atom), and the potential in the interstitial region retains a general form, typically being expanded over plane waves. The method of solving the KS equations adopted in WIEN2k is a successor of augmented plane waves method by [START_REF] Slater | An augmented plane wave method for the periodic potential problem[END_REF], where a composite form of one-electron wavefunctions was applied: numerical radial functions inside the muffin-tin spheres, matched onto plane waves in the interstitial. Following a sequence of modifications, explained in the Introduction to the WIEN2k User's guide (that cites, in its turn, the essential sources), nowadays a combined basis set is applied in the WIEN2k method. It includes numerically tabulated radial functions (centered at each participating atom, for different orbital moment values l) that are matched onto the plane waves at the spheres, plus a number of "local orbitals" which are "relieved" of matching conditions. 1 The number of basis fuctions amounts to "large"

Practical solution of Kohn -Sham equations

(∼thousands; numbered by planewaves with their corresponding cutoff, see the above discussion on VASP) plus "small" (just a few "local orbitals" per atom) hence, generally, "large" as in planewave methods (see the details of the present calculations below in Subsec. 1.3.4). It should be noted, however, that these basis functions are (i) more sophisticated than just plane waves and able to describe strongly fluctuating shape of electron wavefunctions within the atoms, an important improvement over "just plane waves" methods; (ii) adjustable (via their internal numerical part) from one iteration to the other and hence more efficient than "fixed basis" methods. De facto, WIEN2k is able to provide an ultimate numerical precision in practical solving of KS equations, within given approximations adopted at the level of DFT. Therefore this method is often served as an ultimate benchmark for other computational codes; it is used in this quality also in the present work.

Specifics of calculations on periodic systems

It was already briefly mentioned above that the use of planewave basis (or, that of augmented plane waves) implies a periodicity of the system under study: the plane waves are numbered by the vectors of reciprocal lattice G; this presumes that the direct lattice exists, hence the simulation cell must be periodic. However, at this step yet no distinction is done of whether the simulation cell contains an isolated piece of matter (a molecule), or is a part of crystal. The difference between the two cases is the dispersion (of energies, or corresponding wavefunctions) which is non-existent in the first case, but important in the second.

It can be reminded here that dispersion follows from an imposition of the Bornvon Kármán periodic boundary conditions onto the whole crystal. As the consequence, the properties entering the KS equations (1.2) and numbered there by i, a hitherto non specified cumulate index, acquire a more specific numbering in terms of k, the vector within the Brillouin zone (BZ) of the crystal, and the band index ν that runs (upwards) throughout the energy bands, at each k point:

ǫ i → ǫ νk ; φ i (r) → φ νk (r) . (1.23)
As the crystal is assumed infinite, the distribution of k values throughout the BZ is quasi continuous. The summations over electronic states -e.g., the recovering of the electron density according to Eq. (1.1), but also the evaluation of the total energy etc.

-then involve the summation over bands ν and the volume integration over k. For practical purposes, however, only a (more or less) sparse grid of k can be generated, and some kind of interpolation performed. The technics of k-space integration do not depend on the basis set nor the solution Plus Local Orbitals Program", see http://www.wien2k.at/reg_user/textbooks/usersguide.pdf.

Chapter 1 method of KS equations; the most used are the following two approaches, applied with some variations: (i) sampling [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF], whereby a continuous integral of some k-dependent function f is replaced by a weighted sum over its values in N k representative points k i chosen with their corresponding weights w i from within the volume of the BZ, Ω BZ :

BZ f (k)dk → Ω BZ N k N k i f (k i )w i , (1.24)
(ii) the tetrahedron integration, whereby the energy k-dependent properties (energies, matrix elements) are supposed to vary linearly within each tetrahedron delimited by four closely situated k points [START_REF] Blöchl | Improved tetrahedron method for Brillouin-zone integrations[END_REF][START_REF] Jepsen | The electronic structure of h.c.p. ytterbium[END_REF]]. The resulting formulae then depend on the integrand in question, i.e., whether the presumably linear function of k stands in the numerator n(k), in the denominator d(k), or in both, of f (k) = n(k)/d(k), or includes more complex constructions -see, e.g. [START_REF] Molenaar | An extended tetrahedron method for determining Fourier components of Green functions in solids[END_REF].

The both techniques allow to make use of the symmetry relations concerning the integrand f (k) under operations of the space group in question. As a result, the sum (or integration) in (1.24) over the full BZ can be replaced by operations with a reduced set of non-equivalent k points, or over the irredicible part of the BZ. In particular, the weights w i in (1.24) will be affected by taking into account the degeneracies of the (irreducible) k points over which the summation is retained. Whatever the details of implementation, the ultimate accuracy of k-space integration scheme depends on the fineness of the k-grid, that is, essentially (as the homogeneous grids are typically used) on the number of divisions N 1 , N 2 , N 3 along the edges of the reciprocal unit cell. It usually makes sense to choose a "balanced" division according to actual lengths of the reciprocal vectors in question, so that the accuracy of the k-mesh would be governed by a single easily controllable parameter, e.g., the density of k points as discussed in the main text. It is important to make sure that the fineness of the k mesh is sufficient to guarantee the convergence of results to the accuracy relevant for the project in question; such tests are discussed in the following section.

A related but different issue than that of integration over the BZ is that of studying energy dispersions, in view of discussing the placement and composition of valence and conduction bands, etc. For such analysis, band structure plots are usually constructed, by following a quasi continuous k path along the Brillouin zone. Since this serves purely illustrative purposes, a choice of such path is, in principle, arbitrary; however, in order to better emphasize the diversity of the band structure, commonly the path is chosen that connects points of high symmetry of the BZ, including Γ and some points at the BZ edges. The reason for such choice is that the extrema of different 1.3. Practical solution of Kohn -Sham equations bands do usually occur at k points that possess certain symmetry. The established attribution and notation of symmetry points for different types of crystal lattices can be found, e.g., in Chapter 3 "Space Groups" of the book by [START_REF] Bradley | The Mathematical Theory of Symmetry in Solids: Representation Theory for Point Groups and Space Groups[END_REF]. This information is moreover scattered over internet ressources (e.g., http:// lamp.tu-graz.ac.at/ ~hadley/ss2/appendix/symmetrypoints.pdf at the web page of Peter Hadley at the TU Graz), and is organized into tables (without figures, for each given space group) at the Bilbao Crystallographic Server, section "The k-vector types and Brillouin zones of Space Groups": http://www.cryst.ehu.es/cryst/get_kvec. html. A more specific discussion concerning the structures of the present study will be given in Chapter 2. Along with the band energies, other features can be shown at band structure plots, for instance, the selected orbital contributions at a given band and k point. These latter amount to the decomposition of φ νk (r) from (1.23), projecting them onto atom-and orbital character (s, p, d) -related reference functions.

Practical calculation: important files and parameters

In probably any ab initio calculation method, the description of the crystal structure (translation vectors and coordinates of atoms in the unit cell) needs to be explicitly given, accompanied by switches and parameter values typical for the method in question, and according to a particular calculation task: calculating or not the forces on atoms, displacing or not the atoms etc. One should distinguish between the electronic relaxation, that is, finding a self-consistent solution of the KS equations (1.2) for a given fixed geometry, and the ionic relaxation, that consists in finding the equilibrium crystal structure. The latter is achieved by minimizing the total energy (at zero temperature), whereby the lattice vectors and the atom positions in the cell are so modified so as to minimize the forces acting on atoms and the components of the stress tensor. 1 The simplest isolated calculation task is to converge the electronic-structure problem, i.e., to perform the electronic relaxation, for a given (fixed) crystal structure. The selfconsistency of the KS equations implies that the eigenvectors of the KS equations (1.2) give rise, via Eq. (1.1), to exactly the charge density that enters the Coulomb and XC potential of the KS equations. Such convergence is (hopefully) achieved in a sequence of iterations. A usually applied criterion of convergence is such that the charge density deviation from its values at the previous iteration(s), somehow estimated throughout the unit cell, becomes less than some preset value. Setting such a tolerance parameter is reasonably done in both WIEN2k and VASP by default and does not need any special attention. The speed and the stability of the convergence are controlled by the choice of the mixing scheme and mixing parameters. Some technical details will be discussed below.

Chapter 1

WIEN2K

The documentation to the current version of the code can be found at the projet's web site, http://www.wien2k.at/reg_user/textbooks/usersguide.pdf. In a nutshell, the crystal structure of the system under study is defined by the "master input file" case .struct, whereby case stands for the generic name of the calculation task, that must also be the name of the work directory in which numerous working files will reside, and .struct is the fixed extension. A valid example of such file (with large parts omitted and indicated by points) is shown in Fig. 1.1. In the fixed format accepted for this structure file, the line with lattice parameters (a = b = 7.0934 Bohr ≈3.7537 Å, c = 30.458 Bohr ≈16.1177 Å; α = β = 90 • , γ=120 • ) is followed by lines with relative coordinates of constituent atoms: z = ±0.57567 for Ga and ±0.14851 for Se; otherwise (x, y) = ±( 13 , 2 3 ); further Wyckoff positions of this structure are omitted in the example shown. The space group identification (here: Nr 187, P 6m2), local rotation matrices following each nonequivalent species and the symmetry operations (rotation + translation matrices; 12 in the present case) do not need to be correctly indicated NUMBER OF SYMMETRY OPERATIONS -1 1 0 0.00000000 -1 0 0 0.00000000 0 0-1 0.00000000 1 ( ... ) 12 Omitted fragments are indicated by dots in brackets. See text for details.

Practical solution of Kohn -Sham equations

from the beginning since they can be created / verified / updated by corresponding parts of the WIEN2k code. A rather sensitive element of this input file is the muffin-tin (MT) radius, R MT , to be specified for every non-equivalent species (here: 2.08 Bohr for Ga and 2.18 Bohr for Se). Its choice is, at the same time, not critical in the sense that the physically relevant results are not supposed to significantly depend on it (if chosen within some reasonable limits); at the same time its choice is crucial, because once chosen the MT radii cannot be smoothly modified in the calculation. The point is, R MT defines the matching point of the numerically tabulated atom-centered basis function near the atom core with its plane-wave augmentation in the inerstitial region.

A large R MT allows to grasp a larger part of the wave functions / charge density, as spatially fluctuating as they may be, by a more accurate numerics, controlled by a relatively unproblematic and not expensive expansion over spherical harmonics around a given center. A smaller R MT shifts the "load" of accurate numerical description of spatial functions more into the interstitial, demanding a higher planewave cutoff and hence larger matrix sizes to diagonalize. The essential limitation is that the MT radii may not overlap in the course of calculations involving ionic relaxation, displacements of atoms and distortions of the unit cell. Therefore for a calculation project intended to incorporate a dynamical part, the R MT are usually chosen small enough to allow some free space around, to move the atoms. The choice of slightly different set of R MT radii may result in markedly different absolute values (in the program's internal scale) of band energies ǫ νk and the total energy. This is typically not a problem, because the physically relevant properties (e.g., band energies relative to the Fermi level; the defect formation energies or cohesion energies) represent a difference of two "absolute" values which are hardly of interest by themselves. Such absolute energy values, whenever shown in the present work, are only given for reference properties; it is more their variation as function of some parameter (k-mesh density; volume, etc.) that is brought into the corresponding discussion. After some preliminary tests, the R MT radii consistently accepted in the present work are the following: 2.08 Bohr for Ga, 2.27 Bohr for In, 2.18 Bohr for Se, 2.07 Bohr for Cu. They are not advised to be chosen very different for atoms whose sizes do not markedly differ; otherwise the R MT are the user's choice and not subject to any optimisation.

A typical calculation by the WIEN2k code is controlled by a shell script that invokes in a sequence the executable files, created by compilation of (mostly) Fortran95 program files. A single step of electronic relaxation is realized by a sequence of calls to (in principle, stand-alone) programs entitled lapw0 (that constructs the Coulomb and XC potential from the available charge density), lapw1 (diagonalizing of KS matrices), lapw2 (that recovers the valence-electrons charge density from the eigenvectors), lcore (recalculating the core states in each atom), and mixer (that admixes the just obtained charge density to the "old" one, to prepare the input for the following iteration. The Chapter 1 exchange of data between the codes is organized via a system of files which are created or read at corresponding steps. Special shell scripts are provided for preparing necessary files to start a new project, to extract "useful" information (band structures, density of states, charge densities etc.) for inspection. Such organization of calculation is robust, because each step can be separately inspected and, in case of problem, the concerned data files repaired. On the other hand, the organisation of involved calculation tasks might appear cumbersome.

For instance, the full structure optimisation of a given compound in WIEN2k is not easy. A shell script min lapw exists which, for fixed lattice parameters, displaces (retaining the crystal symmetry) the atoms according to forces acting on them. That means that the coordinates in the "active" .struct file are modified (overwritten), till all the forces get smaller than the certain tolerance value. However, the modification of the lattice parameters is not so automatised and demands scanning more or less by hand. In practice, a sequence of trial .struct files can be automatically generated for a given structure, corresponding to a scan through a certain number of lattice parameter values, combinations of the latter -for instance, the c/a and unit cell volume ∼a 2 c for tetragonal or hexagonal structures (as those concerned by the present work). The calculated total energy values through the set of trial structures can be then fitted to a polynomial, and the equilibrium lattice parameters estimated. If these modified lattice parameters happen to be far from those previously assumed, another minimisation of the forces on atoms might be due. The equilibrium structure thus emerges from a sequence of consecutive [lattice parameters optimisation] -[atom coordinates adjustment] steps.

Among the parameters controlling a WIEN2k calculation, major part of which are "safe" in the sense of numerical accuracy at their provided default settings, there are some whose effect on the accuracy of results (and hence on quantitative and qualitative conclusions thereof) might be considerable. The RKMAX parameter is a (unitless) product of the plane-wave cutoff K max by the smallest MT sphere radius in the system under study; its default value of 7.0 is quite acceptable for producing accurate band dispersions of densities of states; however, a good convergence of total-energy results demands a much higher cutoff. We used the value RKMAX=9.0 in all WIEN2k calculations discussed. Another sensitive parameter is the density of the k-mesh for integrations over the Brillouin zone. The mesh practically used (i.e., the number of divisions along the three vectors of the reciprocal cell) is generated after prompting the user for a rough estimate of the desired number of k-points during the initialisation of a WIEN2k calculation. The initialisation sequence, run by the shell script init lapw, invokes the code kgen that takes care of setting up the mesh. In sensitive cases, the convergence of physically relevant total energy differences (e.g., between different crystallographic phases, like e.g. for the case of binary GaSe / InSe structures discussed in the present work) may be very slow. Since this analysis equally concerns the utilisation of the VASP code, the corresponding discussion is postponed till after the brief discussion of the both codes.

Practical solution of Kohn -Sham equations

VASP

Differently from WIEN2k, VASP is the "all-purpose" single code capable to manage electronic and ionic relaxations, the latter to be done with and without distorting the lattice along with modifying the atoms' internal coordinates. The VASP team maintains the project site http://www.vasp.at with links to the exhaustive manual and the "VASP wiki". Among the "standard" fixed-name input files of the VASP code, the POSCAR file is that specifying the crystal structure. A "minimalist" example in Fig. 1.2 refers to the ε structure of GaSe, not literally identical to that shown in Fig. 1.1. The lattice vectors, given line by line with their Cartesian components in Ångström, are followed by relative atom coordinates; four Ga atoms followed by four Se ones. Differently than the all-electron WIEN2k code that generates all the necessary atom-related information "from scratch", VASP relies on POTCAR file, which provides the pseudopotentials for all constituent species, typically from the standard VASP database. The two other necessary fixed-name input files needed for VASP are the KPOINTS which explicitly describes the k-mesh sampling (somehow similar to that used in WIEN2k) and the INCAR file accomodating a variety of parameters that specify the calculation scenario. Many of these parameters are good in default setting, but some need attention:

• The PREC tag, that affects certain cutoffs, is set to "Normal" by default (since version VASP.5) whereas the accuracy required in our calculations demanded (at least) the "Accurate" level.

• ENCUT, the planewave cutoff for the basis functions defined according to Eq. (1.22), needs to be (at least) 500 (measured in eV) for the systems and problems of the present Chapter 1 study.

• EDIFF, the convergence criterion for stopping the electronic relaxation, has to be reduced to, at least, 1E-6 (the default being 1E-4 = 10 -4 ).

• In order to allow the ionic relaxation, the switch IBRION=2 needs to be set up that chooses the conjugate-gradient algorithm (otherwise no relaxation is done by default); moreover EDIFFG = -0.03 stipulates, as the criterion of termination, that forces on all atoms ought to be inferior to 0.03 eV/ Å.

• The ISIF tag chooses the ionic relaxation mode; it provides a number of possibilities to allow or not the displacement of atoms (driven by forces) and variations of the unit cell volume and shape (driven by calculated stress tensor). The "interesting" combinations are: ISIF=2 (atoms relaxed in rigid unit cell, a scenario corresponding to running the min lapw script in WIEN2k) and ISIF=4 (relaxation of ions and of the unit cell shape, for fixed cell volume; a convenient mode for constructing energy/volume curves to compare competitive crystal structures). ISIF=3 performs an unconstrained relaxation of all structure parameters.

The crystal symmetry in VASP is by default taken into account; however, it can be explicitly switched off assigning the ISYM tag a value of 0 or -1. Specifically, some symmetry analysis is done, and the atom coordinates and the charge density are symmetrized. The WIEN2k code enforces a more strict use of symmetry than VASP: a detailed space group analysis is done, and all the relevant calculated properties consistently symmetrized. An attempt to switch off the existing symmetry would demand an arduous trickery (describing equivalent atoms as belonging to different species etc.) whereas in VASP, a modification of a single tag would do the job.

We come now to discussion of the k-mesh convergence issue, that is somehow common for WIEN2k and VASP (and for many other electronic structure codes). The accuracy of k-space integration according to Eq. (1.24) by either the tetrahedron method Blöchl et al. [1994]; [START_REF] Jepsen | The electronic structure of h.c.p. ytterbium[END_REF] or by sampling [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF] improves as the k-mesh becomes more fine. The criterion of convergence in the sense of k-mesh "saturation" is that the properties of interest (typically, total energy differences between different states subject to comparison) become stabilized to the meaningful accuracy.

An example of such analysis is given on Fig. 1.3, for the case of converging the total energies of several phases of GaSe, on the basis of WIEN2k calculation. The density of k-mesh is characterized by a single parameter, the (cubic root of) the number of k points within the unit of reciprocal volume. One can see that the absolute total energy values goes on to decrease as the k-points density (Dk) increases. However, what is of physical meaning is not the total energy by itself, but its difference between the phases one wants to compare. The Dk necessary to get reliable trends may vary from case to case. Obviously, the energy difference between the β and ε phases gets stabilized Total energy per 2 f.u. + 17485000 (mRy) (on the energy scale of ≃ 0.1 mRy, relevant for discriminating the two) already from Dk ≃ 25 Å; the δ phase reliably places itself between the latter phases from Dk ≃ 30 Å on. Finally, the γ phase converges worse and arrives at the same slope of ∆E/Dk only from Dk ≃ 55 Å. In fact, the γ phase has rhombohedral unit cell (and Brillouim zone), whereas the other three phases are hexagonal. Moreover the β and ε phases share roughly the same c/a relation while the δ unit cell is twice longer. Consequently the systematic errors of k-space integration are the most efficiently compensated for β vs ε phases, whereas the γ phase has the tetrahedra completely differently arranged in the k-space, hence such systematic error compensation does not occur. The message from this analysis is that the (N 1 ×N 2 ×N 3 ) divisions along the reciprocal lattice vectors need to be taken as (at least) (17×17×3) for β and ε phases, (17×17×1) for the δ phase and (12×12×12) for the γ phase, in order to yield a reliable energy comparison between the three phases. An alternative solution concerning the "problematic" γ phase could have been to treat it in hexagonal setting, that would permit a not so demanding k mesh Chapter 1 (with Dk ≃ 35 Å), at the price of triplicating the number of atoms in the unit cell. However, in terms of calculation time, a multiplication of the number of atoms is much more demanding than an "equivalent" increase of k points. Moreover, at least when using WIEN2k code that is quite strict in exploiting the maximal available symmetry, it costs efforts to cheat the code into imposing a (lower) hexagonal symmetry onto an inherenty rhombohedral system.

( 1 7 × 1 7 × 3 ) ( 1 4 × 1 4 × 2 ) ( 1 0 × 1 0 × 2 ) ( 9 × 9 × 2 ) ( 7 × 7 × 7 ) (6 ×6 ×6 ) ( 8 × 8 × 8 ) (7× 7×1 ) ( 1 7 × 1 7 × 3 ) (17×17×1) ( 1 0 × 1 0 × 1 0 ) ( 1 2 × 1 2 × 1 2 ) β γ δ ε
As a final remark, Fig. 1.3 serves just test purposes and does not reveal the true energy relation between phases. The calculation compared were done for some fixed geometry of each phase, close to (but not exactly at) corresponding equilibrium structures. A more meaningful inspection of energy relations between phases has to take into account the variation of energy with volume, and is discussed in due place.

Chapter 2 Systems to be studied

General introduction

The initial motivation for the work was to study the electronic structure of mixed semiconductors important for photovoltaic applications. Prototype systems in this relation are chalcopyrite-type CuInSe 2 (CISe) and CuGaSe 2 (CGSe). CISe and CGSe crystals, among many other ternary compounds with chalcopyrite structure, have been probably first described and synthetized by [START_REF] Hahn | Untersuchungen über ternäre Chalkogenide. V. über einige ternäre Chalkogenide mit Chalkopyritstruktur[END_REF]. The practical interest for related systems arised in 1970s in the context of optical applications, e.g., the luminiscence of CuGaS 2 [START_REF] Wagner | Green electroluminescence from CdS-CuGaS 2 heterodiodes[END_REF] or CuInS 2 [START_REF] Bridenbaugh | Junction electroluminescence in CuInS 2[END_REF]]. The corresponding complex sulphides, or mixed sulphides-selenides, do also appear in practical works on photovoltaics. As polycrystalline thin films used for photovoltaics, these systems marked an outstanding efficiency, up to 19.9% [START_REF] Ramanathan | Properties of 19•2% efficiency ZnO/CdS/CuInGaSe 2 thin-film solar cells[END_REF][START_REF] Repins | 19•9%-efficient ZnO/CdS/CuInGaSe 2 solar cell with 81•2% fill factor[END_REF]. CISe and CGSe possess "good" band gap for the use in photovoltaics, that can be further on tuned by alloying [START_REF] Gabor | Band-gap engineering in Cu(In,Ga)Se 2 thin films grown from (In,Ga) 2 Se 3 precursors[END_REF]. [START_REF] Murata | Optimum bandgap profile analysis of Cu(In,Ga)Se 2 solar cells with various defect densities by SCAPs. Japan[END_REF] write: "The bandgap of a Cu(In,Ga)Se 2 (CIGSe) absorbing layer is varied from 1.0 to 1.7 eV by changing the composition ratio of gallium, realizing an optimum design for solar cell absorbers." Many other compounds of the same family have been studied, under the angle of possible photovoltaic applications or otherwise. A rich overview of the knowledge corresponding to the end of 1970s in what concerns structures, band gaps and positioning in the structure phase diagrams of several important chalcopyrite-type ternary semiconductors was done by [START_REF] Pamplin | Ternary chalcopyrite compounds[END_REF]. Concerning first-principles calculations, a pioneering work of [START_REF] Jaffe | Electronic structure of the ternary chalcopyrite semiconductors CuAlS 2 , CuGaS 2 , CuInS 2 , CuAlSe 2 , CuGaSe 2 , and CuInSe 2[END_REF] should be mentioned which served as a benchmark (for band structures, charge densities etc.) or a source for inspiration in more than 360 later publications.

The presence of structural defects, either intentional (e.g., chemical doping) or unintentional (in the course of samples preparation or treatment) is another hot topic of study on these materials. Generally, it somehow reduces the photovoltaic efficiency -Chapter 2 see, e.g., [START_REF] Cao | Defects in Cu(In,Ga)Se 2 chalcopyrite semiconductors: A comparative study of material properties, defect states, and photovoltaic performance[END_REF] and [START_REF] Wi | Photovoltaic performance and interface behaviors of Cu(In,Ga)Se 2 solar cells with a sputtered-Zn(O,S) buffer layer by high-temperature annealing[END_REF], just to mention a couple of recent works. Some defects can be "cured" under appropriate growth conditions [START_REF] Mainz | Annihilation of structural defects in chalcogenide absorber films for high-efficiency solar cells[END_REF], and in general these chalcopyrite-type materials have a reputation to be quite robust, in what regards not much degrading their photovoltaic properties, under irradiation. Among first-principles works, that of Zhang et al. [1998] should be singled out as the one that established some hierarchy in the relative "importance" of possible point defects (in CISe), discussed the placement of the defects' energy levels etc. An extension over point defects in CGSe has been done by the same team soon after [START_REF] Wei | Effects of Ga addition to CuInSe 2 on its electronic, structural, and defect properties[END_REF]].

The defect-related issues have been intentionally left out of the scope of this work, and, for chalcopyrite-type systems, only Ga/In isovalent substitution briefly considered. In view of large number of first-principles calculations done in the course of last decades by many researchers on CISe and CGSe, the present motivation was to try the state-ofart accuracy (within the DFT, using WIEN2k) in optimizing the crystal structure and exploring the band structures, "corrected" within the mBJ XC scheme. To study the effect of Ga/In substitution was part of the work objective; this (partial) substitution was simulated in the simplest way, by considering the number of supercells describing fictitious ordered structures with 1:3, 2:2 ans 3:1 composition of In:Ga. Despite the "crudeness" of this approach, it can be somehow excused by an a posteriori observation that the resulting composition-dependent trends are very smooth.

The chalcopyrite-type structure of CISe (or, similar for CGSe) is shown in Fig. 2.1; the Wyckoff positions of participating elements, corresponding to the space group I 42d (Nr 122), are indicated in Table 2.1. The structure can be understood as the zincblendetype one, in which the cation sites are split into two species, occupied, in an ordered way, by two distinct elements, stemming from the two columns of the periodic table. This lowers the symmetry: it is reduced from cubic (in zincblende) to tetragonal, although the c/2 would only slightly differ from a. Moreover the anion sites acquire a free

Site

Coordinates Element (4a) (0 0 0); ( These experimental data will be discussed later, in comparison with our calculation results.

1 2 0 3 4 ) Cu (4b) (0 0 1 2 ); ( 1 2 0 1 4 ) In / Ga (8d) (x
The nearest coordination of atoms which enter the covalent bonding to form the crystal structure of a semiconductor tend to satisfy the "octet rule", that is, to make a closed 8-electron shell for each atom, counting together the electrons donated by / shared with its neighbours. In case of a conventional (4-coordinated) semiconductor, say ZnSe, whether in zincblende or wurtzite structure, the counting would read as follows:

electr

and similar for III-V (e.g., GaAs), IV-IV (e.g., SiGe) or I-VII semiconductors.

In CISe, each Se atom has 2×Cu + 2×In neighbours, hence 2×1 + 2×3 = 8 valence electrons provided by all these cations, this yet to be divided by 4, because each cation shares its valence electrons to four Se neighbours. Adding 6 electrons provided by the Se itself yields 8. For cations, the octet rule is satisfied "on the average", the electronegativities of Cu and In being not identical. The same counting applies to CGSe, or throughout the compositions in the Cu(In,Ga)Se 2 solid solutions.

As was already mentioned, only hypothetical (ordered) mixed chalcopyrite-type sys- ). This way it is seen that the compound is a mixed selenide, both components of which are "electrically balanced": Cu + 2 Se 2-and In 3+ 2 Se 2- 3 . Obviously, these charge indications are purely nominal, as both compounds have a strong degree of covalency, especially the second one. As n increases, the In 2 Se 3 starts to dominate over Cu 2 Se. The binary compound In 2 Se 3 (and, correspondingly, Ga 2 Se 3 ) have been reported to exist in different phases, whereby the hexagonal one seems to be the most stable. Its crystal structure, however, leave place for some ambiguities. The composition of the hexagonal phase and its discussed strucure models are addressed in Sec. 2.6. An interesting observation is that, on varying the Cu/In concentration throughout ternary systems, a transition must occur at some point from chalcopyritetype (that is, tetragonal) to In 2 Se 3 -type hexagonal lattice. It was one of initial ambitions of the present work to try to suggest a model for this structure transformation, and to compare it with the experimental data which exist on the latice type and the symmetry of some intermediate systems, but do not reveal the placement of atoms [START_REF] Bodnar | Crystal structure of CuIn 3 Se 5 and CuIn 5 Se 8 ternary compounds[END_REF][START_REF] Bodnar | Fabrication and photoelectric properties of oxide/CuIn 5 Se 8 heterojunctions[END_REF]. However, this task has not been fulfilled in the course of work, and only the "pure" In 2 Se 3 , Ga 2 Se 3 systems have been analyzed along this line.

The In 2 Se 3 , Ga 2 Se 3 , in their turn, can be "derived" from "prototype" layered structures InSe and GaSe, correspondingly. These systems are quite interesting, and in fact the major part of the present work is devoted to them. They crystallize in hexagonal structure, that have some similarity with wurtzite. In wurtzite, the sequence of ingerchanged cation/anion hexagonal (0001) layers repeats itself along [0001] in the AB-AB-AB... stacking, so that each given atom, say a cation B, is bonded to three anions in A-layer "below" (no dash) and to one A-atom from the layer "above" (indicated by dash). In InSe / GaSe binaries, one can say that the AB layer is isolated, duplicated as a mirror image, and pinned to its mirror counterpart the cation-cation bonds perpendicular to the layers: AB-BA. The resulting hexagonal double layers are further on stacked in different combinations, giving rise to a variety of polytypes. It can be noted that (i) the outer anion layers are always in the eclipsed configuration, that 2.2. Ordered structures based on cation substitutions in chalcopyrite is, mirror-symmetric with respect to the double-layer mid-plane; (ii) there is no anionanion covalent bond between the outer atoms of two adjacent double layers. Each anion is only bonded to three neighboring cations, whereas each cation is connected to three anions and its opposite (symmetric) cation. Such arrangement helps to formally satisfy the octet rule as follows. In case of a double-layer III-VI semiconductor, 3 own valence electrons per cationic (say Ga) site plus X (yet unknown overlap with the neighboring Ga atom) plus 3 times Y (overlap with each of three neighboring Se atoms) makes 8 to satisfy the octet rule; on the Se site, 6 own electrons are added 3 times Y (overlap to neighboring Ga) to yield the octet rule, hence

3 + X + 3Y = 8 ; 6 + 3Y = 8 ⇒ X = 3 ; Y = 2/3 . (2.2)
In the "fully ionic" counting, the Se atom strips 3×2/3 from its Ga neighbors to fill the own full shell, that leaves 3 -3×2/3 = 1 electron on each Ga atom, to be collectivised on the Ga-Ga bond. We'll see later on how this formal counting corresponds to the bands occupation.

There is, at most, very weak, Se-Se covalent bonding between the adjacent double layers. It is presumed, somehow simplifying, that the layers are kept together essentially by the vdW interactions. In reality the Se-Se bonding seems to be stronger than might be expected from just the vdW interaction. The space between the adjacent double layers is referred to as the vdW gap.

The possible arrangements of double layers is subject to general constraint of hexagonal (or trigonal) symmetry and a couple of simple rules specified in Sec. 2.4. The variety of double-layer stockpilings results in a number of polytypes for a given binary compound. The polytypes relevant for the present work are further on explained in detail, and the works previously done on the binary systems briefly discussed.

Ordered structures based on cation substitutions in chalcopyrite

As the primitive cell of the chalcopyrite structure does already contain two identical In (or Ga) atoms (see Table 2.1), making them different is probably the easiest way to create a In:Ga=1:1 superstructure (Fig. 2.2, left panel). The structure in question will be in fact that of kesterite, named after the Cu 2 (Zn,Fe)SnS 4 mineral but common in new photovoltaic materials -see, e.g., [START_REF] Siebentritt | Kesterites -a challenging material for solar cells[END_REF]. Note that in "genuine" kesterite materials, the cations in the 2,3 positions of the chemical formula typically have split valences, i.e., Zn II + Sn IV rather than In III + Ga III in our formal construction.

Chapter 2

The kesterite structure presumes strict alternance of Cu-Ga and Cu-In (001) planes. However, the (112) planes, the most "densely packed" ones in the chalcopyrite structure [equivalent to (111) in zincblende], are evenly populated with In and Ga atoms. Each Se anion is now within a tetrahedron composed of two Cu atoms, In and Ga, so that no more symmetry constraints apply to keep the Se atom at any symmetric position; all its three internal coordinates are free. The space group of this superstructure is I The suggested 3:1 structure is somehow too artificial in the sense that it imposes a strict (3 Cu-In layers)/(one Cu-Ga layer) ordering, hardly very probable in real materials which would likely favour a more uniform distribution of In versus Ga ions over different lattice planes. However, we retain this simplified model for the sake of easy comparison with pure CISe / CGSe chalcopyrites and with the 1:1 structure.

A consequence of having a simple tetragonal (st) cell for the study of 3:1 structures, in place of a body-centered tetragonal (bct) in case of pure chalcopyrites and a 1:1 substitution, is that the BZ is cut (horizontally and vertically) along the medians of its eight delimiting hexagons, reducing its volume by factor of two. The construction of the corresponding BZs is explained in Appendix 1 (Sec. 2.7), and the resulting shapes, along with (identical) k-path selected for further discussion of band dispersions, are shown in Fig. 2.3. One notes notably that the Z -Γ path in the body-centered structure maps onto the Γ -Z -Γ of the simple tetragonal one, that would yield mirror-symmetric dispersion (half of this path becoming redundant). Another observation is that there must be twice more bands in the (two times smaller) BZ of the bct structure.

Previous works on pure CuInSe 2 and CuGaSe 2 chalcopyrites

The amount of experimental works addressing the optical properties, or otherwise (indirectly) the electronic structure, is tremendous and will not be covered here in detail. Among the theory works that fall more or less in the context of first-principles calculations, mostly DFT, one can single out, as one of the pioneering works, the bandstructure calculation on CuGaSe 2 done by Šipr et al.

[1997], albeit for "simplified" tetragonal lattice. The band dispersions were not specified there; the discussion immediately proceeded towards the interpretation of the X-ray absorption spectra, obtained in the experiment. The "explicit" electronic structures of a number of chalcopyrites, including CISe and CGSe, have been calculated (for fixed experimental lattice parameters) by [START_REF] Jaffe | Electronic structure of the ternary chalcopyrite semiconductors CuAlS 2 , CuGaS 2 , CuInS 2 , CuAlSe 2 , CuGaSe 2 , and CuInSe 2[END_REF]. This study was followed by a detailed discussion, by the same authors, of details of band dispersion and the impact of different XC schemes [START_REF] Jaffe | Theory of the band-gap anomaly in ABC 2 chalcopyrite semiconductors[END_REF].

The next important step concerns the optimization of crystal structure on the basis of DFT calculations. As it seems, the first attempts concerning the systems under discussion have been undertaken by [START_REF] Kılıc | n-type doping and passivation of CuInSe 2 and CuGaSe 2 by hydrogen[END_REF], who used the VASP code Chapter 2 (which experienced important improvements since then), and by Belhadj et al. [2004] who used WIEN97, a predecessor of the WIEN2k code. Both calculations were done with the LDA; predictably, the lattice parameters came out somehow underestimated as compared with experimental data. These results are summarized in Table 2.2. For experimental values, only those cited by [START_REF] Kılıc | n-type doping and passivation of CuInSe 2 and CuGaSe 2 by hydrogen[END_REF] are included in the table. Further experimental results are cited by Belhadj et al. [2004], and a more exhaustive selection can be found in Chapter 3.

In recent times, the implementation of hybrid XC functionals in the calculation codes stimulated a number of new studies, aimed at different problems (related to doping etc.), which included however the results on pristine CISe and CGSe as intermediate benchmarks. Thus, at least three recent works, namely Paier et al. 2.3. A slight but noticeable scattering of the numbers obtained with effectively identical method may give an idea of practical stability / reliability of results that can be considered as "state of art" at the DFT level.

Among "state of art" works done not with the VASP code, one can single out that of [START_REF] Liborio | Atomic structure of the (001) surface of CuGaSe 2[END_REF], dedicated to the study of (001) surface of CGSe, related reconstructions etc. using the CRYSTAL code (an all-electron one, with Gaussian-type basis functions). The results for perfect CGSe bulk can be compared with those in Tab. 2.3 because they are otherwise also done for two XC potentials: (i) the GGA-PBE and (ii) B3LYP, a (different than HSE) hybrid functional. According to this calculation, A considerable "scattering" of first-principles results concerning the lattice parameters and band gaps of pristine CISe and CGSe justifies an attempt undertaken in the present work to probe the XC flavours implemented in all-electron WIEN2k code in what regards their ability to yield reliable structure with PBEsol prescription for the GGA, and to obtain "reasonable" band structures with the mBJ XC functional, more simple in practical calculations than the hybrid functionals otherwise used.

Crystal structures of [In,Ga]Se polytypes

Cation and anion planes are hexagonal, whereby cations are sitting on top of their neighboring cations, and anions are "eclipsed", i.e., along [0001] from the anions of the opposite surface of a given double layer. Obviously, cations are laterally shifted with respect to anions (so as to be at equal distance to the three of them); moreover, the anions are laterally shifted with respect to their counterparts (across the vdW gap) of the adjacent double layer. Remembering that the hexagonal stacking allows three distinct positions, say A/a = (0 0), B/b = ( 2 -Ab-Ba-Ab-Ba-Ab-Baor for the zincblende:

-Ab-Bc-Ca-Ab-Bc-Ca-.

For one of possible II-VI double-layer structures, it may look as follows:

• • • bA-Ab • • • aC-Ca • • • cA-Ac • • • bA-Ab • • •
(dots stand for interactions across the vdW gap), or, in a compact notation:

[Ab Ca Ac] ,
suppressing the redundant information about the mirror symmetry within each double layer, and bracketing the unit of repetition. Only two formal rules have to be respected: (i) cation and anion letters within each given layer must be different; (ii) the anion letters of two adjacent double layers must be different. Apart from these, all combinations are in principle possible, so that, allowing longer repetition chains, one can generate infinite number of polytypes. In practice, only few of them have been identified and seem relevant for the characterization of real systems. The description of different polytypes is scattered over literature. The most systematic source seems to be a work by [START_REF] Kuhn | Crystal structure and interatomic distances in GaSe[END_REF], whose Table 1, with some additions, is reproduced below as Table 2.4. The paper deals with GaSe, so that the internal coordinates may turn out slightly different for InSe.

There is a misprint in the original Table of [START_REF] Kuhn | Crystal structure and interatomic distances in GaSe[END_REF]: the Wyckoff positions for the ε phase enlists Ga in 2(i) position, in place of 2(h). With this error, the ε phase would become a shifted equivalent of β. The Fig. 1 of [START_REF] Kuhn | Crystal structure and interatomic distances in GaSe[END_REF] depicting the structures is correct.

The stacking sequence in Table 1 of [START_REF] Kuhn | Crystal structure and interatomic distances in GaSe[END_REF] is given in terms of just one symbol (letter) per double layer, e.g., the placement of the cation sites, without specifying how the anions are rotated. In this notation, both β and ε structures are coded "ABAB• • • " and hence indistinguishable. 1 We adjust the stacking sequence using the notation explained above; with this it becomes apparent that the β phase reveals "holes" when is looked at along [0001]: one of the three hexagonal (a, b, c) sites remains unoccuopied by neither cations nor anions in none of the planes. For comparison, the differing by the relative rotation (or, translation) of the double layers. The γ phase with triple double-layer stacking has a "virtue" of being a rhombohedral one, hence possessing the primitive cell of just four atoms.

The δ phase is singled out of (many) more-than-three double-layer sequences for the simple reason that it has been included in the crystal structure data by Kuhn et al. repeats as every second one) is different by symmetry from B and C, which are mutually congruent. Neither the placement of the inner cation double plane between the outer anion planes of each double layer, nor the positioning of the B/C cation layers between the adjacent anion "a" layers are fixed by symmetry. This lack of symmetry constraints manifests itself in some distortion (warping along z) of layers, shown in Fig. 2.5 by arrows of largely exaggerated magnitude but correct sense (see Table 2.4). One can note that the warping of layers occurs towards the vacant sites in the adjacent layers above or below. These displacements moreover imply that the z-polarisations within double layers follow in alternating order (those in A layers being opposite to B and C).

The Brillouin zone (BZ) of all phases is hexagonal, with labelling of high-symmetry points as shown in Fig. 2.6. The thickness of the BZ diminishes from β/ε to γ to δ, accordingly to how the c parameter of the unit cell increases. We'll see later (Chapter 4) the similarities and the differences in the band dispersions related to different structures. The case of the γ phase deserves a brief discussion in this particular context (of tracing Note that the choice of Cartesian unit vectors kx, ky differs by 90°-rotation from that of Kx , Ky used in Sec. 2.8. of band dispersions, and comparing them across the phases) . Namely, its electronic structure calculation can be arranged in the hexagonal setting ( similarly to the other phases), with three times the primitive cell (hence 12 atoms) per unit cell; this way the calculations have been done using the VASP code. However, this amounts to spoiling the calculation ressources, because the rhombohedral primitive cell of 'Y contains only 4 atoms and, correspondingly, three timè less occupied bands to take into account in the solution of the KS equations. Moreover, the WIEN2k code, which is particularly strict in applying all the symmetry actually present in the system, strongly "resists" using a unit cell larger than the minimal possibly one. A possible solution would be to perform calculations on the rhombohedral cell and further on, for the sake of producing a "comparable" band structure, to fold down the bands dispersion onto the convenient k path selected in the hexagonal BZ. This "exercise" is explained in the Appendix 2 (Sec. 2.8), specifically for the case of "fiat" rhombohedral BZ (with hexagonal caps, Fig. 2. 7) that corresponds to prolate unit cells ( that is the case for multilayer structures) . In a nutshell, the rhombohedral BZ as in Fig. 2. 7 can be horizontally eut into three slices of equal thickness which (after some cutting/pasting at the edges) become identical to hexagonal prisms, the BZ of the hexagonal lattice. In order to pass from rhombohedral to hexagonal band structure, a path chosen across the rhombohedral BZ must be folded down into a single thin prism. Putting it other way around: in order to reconstruct the band structure in the hexagonal setting from a calculation done in the rhombohedral setting, one needs, by applying the symmetry operations, to triplicate the path chosen within a single hexagonal prism onto the three distinct fragments traced somewhere across the rhombohedral BZ, calculate the band energies along each path, and overlap the results. This way, the resulting band structure will conta.in three times more bands than the initial rhombohedral band structure. An example of such construction will be given in Chapter 4, Fig . .  Before passing on to a. review of previous experimental and theoretical studies of properties, we briefly mention some curiosities related to the crystal structure of III- • GaSe, according to Polian et al. [1976], indeed crystallizes in different modifications, although "the existence of the β-structure, which is well ascertained for GaS, is still doubtful in GaSe" (as of 1976; check later publications). Raman spectra of γ, ε and δ-GaSe are given by Polian et al. [1976].

• Electronic structure of β-GaSe and -InSe was calculated in a tight-binding approach by [START_REF] Camara | Electronic structure of the layer compounds GaSe and InSe in a tight-binding approach[END_REF]. This calculation can be characterized as not a very accurate one, however historically important in view of tight-binding parametrization there elaborated, and due to a comparison with ARPES data done in this work. The work cites, for the crystallographic data on InSe, those from Likforman et al. [1975]. However, Likforman did not refined any phase but γ. Camara et al., for their calculations of β-InSe, took Likforman's a and 2 3 × of their c, and cite the latest number as genuine, without saying a further word...

• In a number of publications, it is repeated that "Indium selenide grown by the Bridgman method crystallizes in the γ polytype" -e.g., [START_REF] Manjón | Experimental and theoretical study of band structure of InSe and In 1-x Ga x Se (x < 0.2) under high pressure: Direct to indirect crossovers[END_REF],and modern calculations are overwhelmingly done for this phase.

• However, Larsen et al. [1977] in their experimental work dealing with band structures extracted from ARPES experiments, specify: "The unit cell of InSe contains eight atoms, two In and two Se atoms from two adjacent sandwiches", not giving any further clue as to whether their single-crystal samples were of presumably β, or ε, type. The comparison to theory in the work cited by [START_REF] Larsen | Determination of the valence band structure of InSe by angle-resolved photoemission using synchrotron radiation[END_REF] referred to then avaiable calculation results for GaSe Schlüter et al. [1976].

• The lattice parameters a=4.05 Å, c=16.93 Å, apparently coming back to Semiletov (1958), left a long incognito trace in the literature. Yu et al.

[1999] cite them as those from Larsen et al. [1977], who, in their turn, cite them as stemming from the Wyckoff's book "Crystal Structures" of 1963.

• The confusion, if not the controversy, concerning the structure of InSe seems to persist in later experimental works. In the combined experimental / theoretical work by [START_REF] Yu | Spin resolved photoemission spectroscopy from InSe(0001) using circularly polarized radiation[END_REF], it was indicated that "InSe is a semiconductor with layered crystal structure belonging to the non-symmorphic space group D 4 6h ", that means P 6 3 /mmc describing the β phase. Chapter 2

GaSe and InSe: previous research of properties

Looking back at the history of theory studies of III-VI systems (mostly GaSe and InSe, some papers also mentioning GaS in the same context), we can identify several "waves" of interest, marked by different motivation.

The end of 1970s saw an initial bunch of band structure calculations, mostly parameterdependent ones, which helped to elucidate band dispersions, and, at least preliminarily, to addres the differences brought about by stacking of layers in different polytypes. Schlüter et al. [1976]; [START_REF] Schlüter | Valence-band density of states and chemical bonding for several non-transition-metal layer compounds: SnSe 2 , PbI 2 , BiI 3 , and GaSe[END_REF] calculated band structure and related properties of GaSe using phenomenological pseudopotentials (PP) and plane-waves (PW) basis. [START_REF] Mccanny | The band structures of gallium and indium selenide[END_REF] calculated band structures of GaSe and InSe using a semi-empirical tight-binding (TB) method (a single tetralayer was taken, hence no distinction between polytypes considered). [START_REF] Doni | Electronic properties of the III-VI layer compounds GaS, GaSe and InSe. I: Band structure[END_REF] performed band structure calculations of β-GaSe and -InSe with a PP TB method. Depeursinge et al.

[1978] calculated electronic structure of β-InSe by semiempirical TB method, and of ε-InSe by an empirical PP method. [START_REF] Depeursinge | Polytypism and layer-layer interaction in the III-VI layer semiconductors[END_REF] discussed band dispersion and folding depending on layer stacking in polytypes (β, ε, γ), on an example of GaSe. "Improved" TB calculation results for β-, γ-and ε-GaSe have been reported by [START_REF] Nagel | Tight-binding study of the electronic states in GaSe polytypes[END_REF]. [START_REF] Robertson | Electronic structure of GaSe, GaS, InSe and GaTe[END_REF] calculated (also within a TB method) the band structures and densities of states of β-GaSe, β-InSe and single tetralayers of both materials, discussing at some length the symmetries and the composition of bands. Some of the mentioned works made a comparison of calculated band structures with angle-resolved photoelectron spectra, e.g., of [START_REF] Larsen | Determination of the valence band structure of InSe by angle-resolved photoemission using synchrotron radiation[END_REF] for β-InSe, and Thiry et al. [1977] for β-GaSe. A much later work by [START_REF] Camara | Electronic structure of the layer compounds GaSe and InSe in a tight-binding approach[END_REF], dedicated to the electronic structure of both InSe and GaSe, also follows the pattern of model TB approach, guided by a comparison to photoemission data. The TB parametrization was such that effectively neglected the difference between polytypes.

Since about mid-1990s, the III-VI semiconductors regained interest for calculations, this time for parameter-free ab initio ones done in the framefork of the density-functional theory (DFT), which proven able to provide accurate structure parameters and meaningfully compare results for different polytypes. Among the methods of practically performing such calculations, some rely on ab initio PPs and often (yet not exclusively) use the PW basis sets. Otherwise, all-electron schemes were used, including the WIEN2k code, with different flavours of the XC potentials used. Most of these works included the relaxation of corresponding crystal structures following the total energy minimization. Some works addressed specific 3-dimensional phases, whereas the others explicitly studied the effect of augmenting the number of (tetra)layers on the 2-dimensional band structure. Gomes da Costa et al. [1993] calculated the electronic structure of β-and γ-InSe by ab initio PP-PW method using the LDA, notably taking into account the spin-orbit interaction. Adler et al. [1998] studied ε-GaSe by a PW method using ultra-2.5. GaSe and InSe: previous research of properties soft PP, within the LDA, whereby the structure optimization was followed by an analysis of phonon dispersion. [START_REF] Yu | Spin resolved photoemission spectroscopy from InSe(0001) using circularly polarized radiation[END_REF] studied spin-resolved ARPES from InSe using circularly polarized radiation; comparison was done with calculations by WIEN97 (the predecessor of the WIEN2k package) taking the spin-orbit interaction into account, however optimizing only the internal coordinates at fixed lattice parameters. Zhang et al.

[2006] optimized the lattice parameters of ε-GaSe using the WIEN2k method and the LDA. [START_REF] Rybkovskiy | Size-induced effects in gallium selenide electronic structure: The influence of interlayer interactions[END_REF] studied band structure of GaSe with 1,2,3,4 tetralayers, using the PP-PW method as implemented in the abinit code, and the LDA. Ghalouci et al.

[2013] calculated the equations of state of β-and ε-GaSe in comparison with other phases (typical for "conventional" semiconductors but too high-energetic and hence "irrelevant" for the III-VI systems), using the WIEN2k method in combination with the GGA. [START_REF] Ma | Tunable electronic and dielectric behavior of GaS and GaSe monolayers[END_REF] provided band structure calculations of GaSe (using the VASP code and GGA) as bulk crystal and as a multilayer system (with 1 to 4 tetralayers) and reported the corresponding lattice parameters. [START_REF] Olguín | Ab initio electronic band structure study of III-VI layered semiconductors[END_REF] relaxed the structure of γ-InSe and ε-GaSe, also using the WIEN2k with GGA (with some additional efforts in the study of band gaps, see below). Zhang et al.

[2014] calculated equilibrium structure, elastic and optical properties of ε-GaSe by a PP-PW method within the LDA. Rybkovskiy et al.

[2014] calculated band structures of β-GaSe, InSe, and GaS in dependence on number of tetralayers, using another realisation of the PP-PW formalism (Quantum Espresso) and several flavours of GGA for comparison, with spin-orbit taken into account.

A number of relatively recent works primarily addressed the properties (optical, elastic, vibrational) under hydrostatic pressure, or under stress, often in the context of comparison or combination with experimental studies. Even as the simulation under pressure are not by itself our interest in the present study, we mention such works below for possible references, since we span the range of volumes while constructing the equation of state curves around equilibrium; moreover some of the "under pressure" works contain the useful reference data concerning the ambient conditions. Manjón et al. [2001] simulated, in combination with experiment, an evolution of some relevant transitions in γ-InSe with pressure, whereby the calculations have been done by the Siesta code, that relies on norm-conserving PPs and basis function of localized numerical functions, using specifically the LDA. [START_REF] Ferlat | Ab initio studies of structural and electronic properties of solid indium selenide under pressure[END_REF] applied several ab initio methods to study γ-InSe at 0, 7 and 8 GPa. [START_REF] Rushchanskiȋ | The influence of hydrostatic pressure on the static and dynamic properties of an InSe crystal: A first-principles study[END_REF] optimized the structure of γ-InSe in a PW calculation (using the abinit code), for a range of pressures from 0 to 14 GPa. [START_REF] Errandonea | Crystal symmetry and pressure effects on the valence band structure of γ-InSe and ε-GaSe: Transport measurements and electronic structure calculations[END_REF] calculated electronic structure (bands in the vicinity of the gap) for γ-InSe and ε-GaSe (hexagonal) at 0 and 7 GPa, discussing the comparison with measured transport properties. [START_REF] Schwarz | Effect of pressure on the structural properties and electronic band structure of GaSe[END_REF] studied structural parameters of GaSe under pressure, whereby the electronic structure of ε-GaSe was calculated by WIEN97 (a predecessor of WIEN2k) in combination with GGA. [START_REF] Errandonea | High-pressure electronic structure and phase transitions in monoclinic InSe: X-ray diffraction, Raman spectroscopy, and density functional theory[END_REF] studied high-pressure monoclinic and tetragonal phases of InSe, applying, Chapter 2 in combination with first-principle calculations, the X-ray diffrection and Raman spectroscopy. In both these works, the Siesta code was employed for calculation. [START_REF] Brudnyi | Charge neutrality level and electronic properties of GaSe under pressure[END_REF] calculated lattice parameters and electronic structure of ε-GaSe (using PP-PW formalism as implemented in the abinit code, in combination with the LDA), also under pressure. [START_REF] Zhu | Topological phase transition in layered GaS and GaSe[END_REF] simulated variation of band structure of β-and ε-GaSe at different strains, using WIEN2k and GGA, in the context of strain engineering of band gaps. [START_REF] Kosobutsky | Structural, elastic and electronic properties of GaSe under biaxial and uniaxial compressive stress[END_REF] calculated electronic and elastic properties of ε-GaSe under biaxial and uniaxial stress, using abinit and Quantum ESPRESSO codes with LDA and GGA. Band structure was compared with ARPES experiments of [START_REF] Plucinski | Electronic band structure of GaSe(0001): Angle-resolved photoemission and ab initio theory[END_REF]. [START_REF] Ma | Tunable electronic and dielectric behavior of GaS and GaSe monolayers[END_REF] discussed the variation, depending on strain, of the band gap in bulk and multilayer GaSe, from the results of band structure calculations done with the VASP code using the GGA.

Several ab initio studies addressed the GaSe compound containing point defects (we are not aware of similar calculations for InSe). [START_REF] Rak | Electronic structure of substitutional defects and vacancies in GaSe[END_REF] simulated electronic structure of doped (substitutional impurities and Ga / Se vacancies) GaSe in a supercell approach, using WIEN2k and VASP code, making reference to calculation of the perfect bulk crystal. [START_REF] Ao | Functionalization of a GaSe monolayer by vacancy and chemical element doping[END_REF] calculated (by VASP method, using the GGA) the electronic structure, relaxation pattern, formation energies etc. related to point defects: the Ga vacancy and larger "holes" (missing GaSe 3 and Ga 2 Se 6 fragments) in a single tetralayer, using 4×4 or 7×7 supercells.

As is the case in many studies of semiconductors, the issue of the band gap, routinely underestimated in "conventional" DFT (LDA or GGA) calculations, was subject to special attention in some recent works. In principle, the GW approximation [...] seems to be a "state of the art" in this sense; a practical alternative can be hybrid functionals, available along with some calculation codes, or e.g. (not as accurate but much less demanding in implementation and ressources) the modified Becke -Johnson (mBJ) flavour of the exchange-correlation functional, implemented in the WIEN2k code. [START_REF] Ferlat | Ab initio studies of structural and electronic properties of solid indium selenide under pressure[END_REF] reported GW results for some high-symmetry k-points in the Brillouin zone (BZ) of γ-InSe. [START_REF] Rybkovskiy | Size-induced effects in gallium selenide electronic structure: The influence of interlayer interactions[END_REF] performed GW calculations at the BZ center of GaSe (for increasing number of tetralayers and the bulk). [START_REF] Olguín | Ab initio electronic band structure study of III-VI layered semiconductors[END_REF] compared band structures calculated for γ-InSe and ε-GaSe using the WIEN2k code with GGA and mBJ, and discussed the latter results in the context of available GW calculations. [START_REF] An | Systematic investigation on topological properties of layered GaS and GaSe under strain[END_REF] addressed the issue of band gap, as obtained by different methods for ε-GaSe (and also in β-GaS), offering an overview across other available results and implementations. Original calculations have been done by WIEN2k, using LDA, GGA, mBJ, and an implementation of the GW formalism within an all-electron approach, interfaced with WIEN2k; moreover, some GW calculations were done using the VASP code.

Since about early 2000s, a number of works have been inspired by lowdimensionality, nanoscale tunability or other "special" properties of III-VI semiconduc- To conclude the overview of previous works, we mention some dedicated to lattice dynamics, which usually specify the ground-state structure as a benchmark, even if vibrations as such are not within the scope of the present study. [START_REF] Wieting | Interlayer bonding and the lattice vibrations of β-GaSe[END_REF] studied vibrations in β-GaSe by infrared (IR) and Raman specroscopy, and offered an extensive analysis of modes. [START_REF] Hayek | Splitting and coupling of lattice modes in the layer compounds GaSe, GaS, and GaSe x S 1-x[END_REF] reported the variation of Raman spectra of GaS x Se 1-x with concentration. [START_REF] Kuroda | Lattice-dynamical and photoelastic properties of GaSe under high pressures studied by Raman scattering and electronic susceptibility[END_REF] studied vibrations in GaSe (among other compounds; under pressure) by Raman spectroscopy. [START_REF] Julien | Far-infrared spectra of indium selenide single crystals[END_REF] studied vibrations in γ-InSe from far-IR spectra. [START_REF] Adler | First-principles investigation of the lattice dynamics of ǫ-GaSe[END_REF] calculated (and discussed in comparison with experiments) the phonon dispersions in ε-GaSe, using the density-functional perturbation theory (DFPT) and LDA. [START_REF] Rushchanskiȋ | The influence of hydrostatic pressure on the static and dynamic properties of an InSe crystal: A first-principles study[END_REF] calculated, using a PP-PW method (abinit), the shift of vibration frequencies with pressure in γ-InSe. [START_REF] Allakhverdiev | Lattice vibrations of pure and doped GaSe[END_REF] measured Raman and IR spectra of pure and doped GaSe and offered a comparison with own and a number of previous calculations. [START_REF] Rybkovskiy | Ab Initio electronic band structure calculation of two-dimensional nanoparticles of gallium selenide[END_REF] calculated band dispersions in GaSe (tetra)layers of increasing thickness (1 to 4) by abinit method, using the LDA. [START_REF] Lei | Synthesis and photoresponse of large GaSe atomic layers[END_REF] described the growth of large thin GaSe layers and their characterization by vibration spectroscopy (Raman), depending on thickness.

Crystal structure of In 2 Se 3

As was mentioned above, the In 2 Se 3 phase, according to nominal valences of constituents, is electrically balanced. Therefore, the octet rule can be satisfied with the tetrahedral coordination of cations and anions, which must however include the presence of ordered cation vacances (one per two In atoms) to provide the correct stoichiometry:

In III 2 Se II
3 . The question is, how are the In vacancies distributed. Probably all structure studies report different phases of either hexagonal or rhombohedral type. The study of [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF] summarizes previous reports (see Fig. 1 thereof) and suggests casting of all reliable information into two phases: the "vacancies ordered in screw form" (VOSF) and the "layer structure phase" (LSP). Both are basically wurtzite, with different ordering of vacancies. Chapter 2

Vacancies ordered in screw form

The VOSF phase corresponds to √ 3 × √ 3 × 3 multiplication of the wurtzite cell and hence has 9×4 -6 = 30 atoms. (We remember that the wurtzite primitive cell has two cations in the A and B hexagonal positions, with an anion on top of each, ideally at z = 3/8). Fig. 2.8 depicts the hexagonal plane with wurtzite primitive cell spanned by the a w 1 and a w 2 vectors, and lattice vectors of the supercell: a 1 = a w 1 -a w 2 , a 2 = 2a w 2 +a w 1 . The third vector is tripled and transverses 6 (say, cationic) planes. The composition of the latter is shown in Fig. 2.9.

The lattice parameters identified by [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF] are: a=7.14 Å; c=19.38 Å. The placement of vacancies corresponds to Fig. 2 of [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF]. It is seen that, coming from plane to plane, the vacancies make a screw around the ( 2 3 1 3 ). Atom position in the VOSF phase, read in directly from Fig. 2.9, are summarized in Table 2. 5.

In fact, the atoms have to be displaced from these positions due to reduced symmetry of their environment. [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF]] specify the space group as P 6 1 (Nr. 169) or P 6 5 (Nr. 170), which differ in the direction of screw (left-hand or right-hand). The Wyckoff positions are of type (6a) only, as follows:

(x y z), (x-y x z + 1 6 ), (ȳ x-y z + 1 3 ), (x ȳ z + 1 2 ), (x+y x z + 2 3 ), (y x+y z + 5 6 ) for P 6 1 ; (x y z), (y x+y z + 1 6 , (x+y x z + 1 3 ), (x ȳ z + 1 2 ), (ȳ x-y z + 2 3 ), (x-y x z + 5 6 ) for P 6 5 . The first choice seems, indeed, to match the coordinates in Table 2.5, if the origin of coordinates in the (x, y) plane is shifted into the projection of the screw axis, ( 2 3 1 3 ). We easily see that all atomic positions in Fig. 2.9 are indeed general (not symmetric) ones. For example, an atom A at ( 1 3 1 3 0) in its plane is equally surrounded by three vacancies, so its projection seems to be fixed by symmetry. However, the placement of its neighbours in the upper plane (z = 1/6) impose only a mirror symmetry, i.e., the Chapter 2

2.8 Appendix 2. Relation between rhombohedral and hexagonal settings for the Brillouin zone and the band dispersions for the γ phase.

We express the reciprocal lattice vectors a * , b * , c * in terms of so far undefined A, C parameters which will be specified further on. The choice of reciprocal lattice vectors in the Figure is the following:

a * = A 2 , - A √ 3 2 , C , b * = A 2 , A √ 3 2 , C , c * = -A, 0, C . (2.3)
The resulting BZ is shown in Fig. 2.12 in two projections. 1 The positions of some 1 We note that at C/A≥1/ √ 2 the upper hexagon disappears; the case symmetric points can be read from the figure. In the cartesian coordinates (along the axes K X , K Y , K Z in the figure) they are:

C/A = 1/ √ 2 corresponds to
1 : 0, 0, 3C 2 2 : A 2 , 0, C 3 : A 2 -C 2 A , 0, 3C 2 4 : A 2 + C 2 2 , 0, C 2 5 : -A 2 , 0, C 2 6 : -A 2 + C 2 A , 0, 3C 2 7 : -A 2 , A 2 √ 3 + 2C 2 A √ 3 , C 2 8 : -A 2 + C 2 A , A 2 √ 3 -2C 2 A √ 3 , C 2 9 : -A 2 -C 2 2A , 0, 0 (2.4)
The transformation to relative coordinates k x , k y , k z (in units of reciprocal lattice vectors):

   k x k y k z    =    1 3A -1 A √ 3 1 3C 1 3A 1 A √ 3 1 3C -2 3A 0 1 3C       K X K Y K Z    ;
This matrix contains linewise the translation vectors in the direct space, whence the identification of A, C should be straightforward. In fact they are: A = 2/(a H √ 3), C = 1/c H . The resulting (k x , k y , k z ) coordinates of symmetric points are:

1 : 1 2 , 1 2 , 1 2 2 : 1 2 , 1 2 , 0 3 : 2 3 -C 2 3A 2 , 2 3 -C 2 3A 2 , 1 6 -2C 2
3A 2 4 :

1 3 + C 2 3A 2 , 1 3 + C 2 3A 2 , -1 6 -2C 2 3A 2
5 : 0, 0, 1 2 6 :

1 3 + C 2 3A 2 , 1 3 + C 2 3A 2 , 5 6 -2C 2 3A 2 7 : -1 6 -2C 2 3A 2 , 1 6 + 2C 2 3A 2 , 1 2 8 : 1 6 + 2C 2 3A 2 , 1 2 , 5 6 -2C 2 3A 2 9 : -1 6 -C 2 6A 2 , -1 6 -C 2 6A 2 , 1 3 + C 2 3A 2
(2.5)

Discuss now how the rhombohedral BZ is related to a hexagonal one. The direct (a, b, c) and reciprocal (a * , b * , c * ) vectors in rhombohedral (R) and hexagonal (H) settings, in terms of hexagonal lattice parameters a and c, are as follows:

( a H b H c H ) = ( e 1 e 2 e 3 )    a √ 3 2 0 0 -a 2 a 0 0 0 c   
(lattice vectors columnwise);

(2.6) cubic primitive lattice.

Chapter 2 (2.7)

   a * H b * H c * H    =    2 a √ 3 0 0 1 a √ 3 1 a 0 0 0 1 c      
( a R b R c R ) = ( e 1 e 2 e 3 )    a √ 3 6 a √ 3 6 -a √ 3 3 -a 2 a 2 0 c 3 c 3 c 3    (lattice vectors columnwise); (2.8)    a * R b * R c * R    =    1 a √ 3 -1 a 1 c 1 a √ 3 1 a 1 c -2 a √ 3 0 1 c       e 1
e 2 e 3    (reciprocal lattice vectors linewise). (2.9)

We note that the projections of the reciprocal lattice sites onto the K X , K Y plane is about the same for R and H lattices (see Figure ). The R-BZ (shown in thin lines in the figure) is "simplified" into a prism (thick lines), and is flattened by factor of 3. The coordinates of high-symmetry points, indicated as for the standard hexagonal setting, are as follows: 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 6 NUMBER OF SYMMETRY OPERATIONS ...

in K X , K Y , K Z in a * H , b * H , c * H in a * R , b * R , c * R Γ 0 0 0 0 0 0 0 0 0 M A 2 0 0 1 2 0 0 1 6 1 6 -1 3 M ′ A 4 A √ 3 4 0 0 1 2 0 -1 6 1 3 -1 3 K A 2 A 2 √ 3 0 1 3 1 3 0 0 1 3 -1 3 A 0 0 C 2 0 0 1 2 1 6 1 6 1 6 L A 2 0 C 2 1 2 0 1 2 1 3 1 3 -1 6 H A 2 A 2 √ 3 C 2 1 3 1 3 1 2 1 6 1 2 -

Chapter 3

Results for Cu(In,Ga)Se 2 systems

Introduction

As becomes obvious from Sec. 2.3, quite a number of first-principles calculations has been done on pristine CISe and CGSe. Our new contribution is (i) the use of PBEsol XC potential, proven to be quite able for quite accurate predictions of equilibrium structures in many compounds with different levels of ionicity / covalence; (ii) the use of mBJ "meta-GGA", that was shown capable, at least, to "improve" the calculated band gaps for a number of semiconductors and insulators; (iii) to trace the above properties throughout the range of systems with partial In/Ga substitution (even if only several small high-symmetry supercell have been taken into account). In the following, we briefly discuss the results concerning the crystal structures (the PBEsol part) and the band structures (the mBJ part) for Cu 2 In 2 Se 4 , Cu 4 In 3 GaSe 8 , Cu 2 InGaSe 4 , Cu 4 InGa 3 Se 8 and Cu 2 Ga 2 Se 4 supercells in comparison.

Crystal structures optimized with GGA

The first-principles calculations on chalcopyrite-type systems were done with WIEN2k method, using several prescriptions for the XC potential. The structure optimization was performed in a sequence of consecutive optimizations of volume (with fixed c/a) and the c/a value (with fixed volume), as the full automatic optimization of lattice parameters is not implemented in the WIEN2k code. On the contrary, the internal coordinates, for lattice parameters fixed, can be searched for by minimizing the forces of atoms. To this end, a script is provided along with the WIEN2k code.

Table 2.2 and Fig. 3.1 imply that the LDA results underestimate the lattice parameters of both CISe and CGSe by 0.5 -1.5%; in order to assess the previous results in this context, we should indicate that there is, in fact, a considerable scattering of Chapter 3 data concerning each of these pristine materials (due to samples preparation, samples quality, characterization methods, etc.) An idea of such scattering can be given by Fig. 3.1, that "visualizes" the data of Table 2.2 in combination with the selection of experimental lattice parameters, Table 3.1.

As concerns the calculation procedure, it can be reminded here that the structure optimization done by WIEN2k is organized step by step in a practically "by hand" way, because the optimization of the lattice parameters is separated from the optimization of internal coordinates. The first is done by fitting the total energy values from certain number of trial calculations with different a, c values. Once the lattice parameters which minimize the energy are identified, the script mini lapw is started that consequently shifts the atoms (within the fixed unit cell) in a loop over electronic structure calculations, until the forces on all atoms become smaller than the specified tolerance value. After which, a refinement of a, c values in another set of scanning the total energies might be needed etc. This makes the procedure cumbersome and prone to human errors... As has been mentioned in Chapter 2, there is only one free internal coordinate (x for anions) in the chalcopyrite structure but more free coordinates in the supercells that describe intermediate concentrations. In fact, the anion is not fixed by symmetry in the center, or along the symmetry axis, of a tetrahedron formed by different cations. Chapter 3

Crystal structures optimized with GGA

Moreover, in the supercells representing the 3:1 (or 1:3) concentration of In:Ga, the "majority" cation can relax along the z direction. The Table 3.2 summarizes the structure information concerning the relaxed end-composition and intermediate-composition structures.

The variation of lattice parameters is also depicted in the left panel of Fig. 3.2. In the right panel of the same figure, a selection of experimental data done by [START_REF] Schnohr | Compound semiconductor alloys: From atomic-scale structure to bandgap bowing[END_REF] is shown. Concerning these trends, one can make the following observations:

• The slope is very smooth and linear with concentration (hence following the Vegard's law, as is common for many mixed semiconductors).

• The tetragonal distortion is very small (c/(2a)≈1 within 1%); however, with concentration a shift occurs from slightly prolate CISe to slightly oblate CGSe.

• The agreement with experimental lattice constants is very good for CISe (note that the old c value reported by [START_REF] Hahn | Untersuchungen über ternäre Chalkogenide. V. über einige ternäre Chalkogenide mit Chalkopyritstruktur[END_REF] is at variance with the later measurements, e.g., those shown in the right panel of Fig. 3.2). However, the lattice relaxation of CGSe yields a bit too short a and too long c (as was already seen in Fig. 3.1). This tendency for "overbonding" for gallium will be seen again in the next chapters concerning binary compounds.

• The crossing point a = c/2 is shifted in the calculation towards ∼75% of Ga whereas Table 3.2: Crystal structure parameters from the relaxation of Cu(In,Ga)Se 2 supercells using the WIEN2k code and PBEsol-GGA. Only "non-trivial" internal coordinates are indicated. in the experiment it is placed at ∼25% of Ga.

Space

The variation of the internal coordinates can be otherwise characterized by the modification, and corresponding scattering, of bond length, as shown in Fig. 3.3. Whereas the lattice constants change by ∼3% throughout the concentration range, each species of the bond length follows this general trend but within much reduced range, within ∼0.7%. Such relative invariance of the individual (element-specific) bond lengths in semiconductor alloys is a since long established trend, known from early EXAFS experiments; a number of examples can be found in a review paper by Claudia Schnohr [2015]. In this work, the bond lengths in the CISe-CGSe mixed systems have been measured; the results are directly compared with our calculations in Fig. 3.3. One should pay attention that the ordinate axis in the left panel of Fig. 3.3 is split into two pieces, so that the In-Se bonds are in fact much longer (by ∼0.2 Å) than the other two bond species which differ by merely ∼0.05 Å (according to calculation) or by ∼0.01 Å (according to experiments). Putting it differently, the Se atoms within the tetrahedra composed of Ga and Cu atoms tend to go more off-center according to calculation than as is derived from EXAFS experiments. Apart from this observation, the agreement between the experiment and theory is quite good. However, one more point might deserve a more careful analysis in the later works. The bond length values for every given concentration are characterized, on the basis of EXAFS data, by a single value. [START_REF] Schnohr | Compound semiconductor alloys: From atomic-scale structure to bandgap bowing[END_REF] specifies (in the caption to her original Fig. 41) that the uncertainty of the values is ±0.002 Å. The problem is, however, that in mixed systems some scattering Chapter 3 of the bond length values occurs simply because, as the local symmetry is lowered, the Se atoms tend to move off-center in all senses, placing themselves closer to Ga than to In. In very ordered supercells like those used in our calculation, this gives rise to a splitting of individual bond species into three distinct In-Se ones in the Cu 4 In 3 GaSe 8 supercell (or, symmetrically, three Ga-Se bond length values in Cu 4 InGa 3 Se 8 ), and four different Cu-Se bonds in each of these cases. In "real" solid solution (say simulated by a big number of large supercells) this scattering will be more smooth, but the general span of the bond length values is not likely to be much reduced. From the left panel of Fig. 3.3 we can estimate the half-width of this scattering as ∼0.1 Å. The analysis of the EXAFS data, apparently, cannot resolve the "fine structure" in the bond lengths distribution and extracts simply the mean value. It may be that the inspection of vibration spectra, which are known to be quite sensitive to bond lengths distribution, would help to provide additional information. Such works on mixed semiconductors (not chalcopyrites so far) are underway in the host laboratory; a more specific discussion is beyond the scope of the present work.

Band structures of pristine systems calculated

with GGA and mBJ

PBEsol-GGA

The band structure of CGSe and CISe is well known from previous theory works, some of which are mentioned in Chapter 2. The study by Soni et al.

[2011] might be particularly useful, because it discusses in detail the band structure and partial densities of states of both systems, as obtained with two calculation methods, one of which is WIEN2k, as in our case, the only slight difference being the XC flavour (PBE in the work cited vs PBEsol in the present work). To set a frame for the following discussion, Fig. 3.4 shows the band structures for CGSe and CISe according to PBEsol. The narrow band near -15 eV is due to Ga3d or In4d states; the isolated band just above it, around -14• • • -13 eV -the Se4s band. The following isolated band near -7 eV is due to Ga4s or In5s states. Immediately above it, separated by a narrow gap, starts the "main" valence band, formed by Cu3d, Se4p and Ga4p or In5p states, within which one finds yet another pronounced gap (at ∼ -2.5 eV), probably due to splitting into bonding and antibonding states. One can note that this latter gap does not come out in all calculations; for instance, it does not appear (neither for CGSe nor for CISe) in the CRYSTAL03 calculation (using the Gaussian basis set) of [START_REF] Soni | Electronic and optical modeling of solar cell compounds CuGaSe 2 and CuInSe 2[END_REF]. At the same time, the WIEN2k results from the same work are practically identical with ours, in what concerns the band structures. The counting of bands, including semicore states, proceeds from counting contribut-3.3. Band structures of pristine systems calculated with GGA and mBJ ing electrons (excluding core shells) as follows: Cu 3p 6 3d 9 4s 2 → 17; In 4d 10 5s 2 5p 1 → 13; Ga 3d 10 4s 2 4p 1 → 13; Se 3d 10 4s 2 4p 4 → 16, hence 17 + 13 + 16×2 = 62 valence-band electrons per formula unit, or 62×2 = 124 electrons per unit cell, or 124/2 = 62 bands. We find this counting -the last occupied band is #62 and the first vacant one #63 -also in Fig. 2 and 3 of Soni et al. [2011]. The Cu3p and Se3d states are deep and not included in the Fig. 3.4. The "standalone" semicore states (up to Ga 4s / In 5s) include 42 bands, so that there rest 20 for the "main valence band" (Cu3d 4s + Se4p + cation p).

Comparing band structures of CGSe and CISe, one notes (i) slightly deeper placement and smaller dispersion of the Ga3d band with respect to the (less localized) In4d band, (ii) a slightly more shallow placement and smaller width of the Se4s band in CISe than in CGSe (probably, primarily a consequence of larger lattice constant in the first compound), and (iii) a much more shallow placement of the In5s band compared to Ga4s. Apart from this, there is a striking similarity between CGSe and CISe in what Chapter 3 concerns dispersions throughout the "main valence band", with one marked exception: at the valence band top in Γ, three upper bands converge and become nearly degenerate in CGSe, whereas in CISe the band #60 is shifted down by ∼0.4 eV. Exactly the same behavior is seen in Fig. 2,3 of Soni et al.

[2011] without being discussed; strangely, in their CRYSTAL02 calculation the three topmost valence bands nicely converge towards Γ for both CGSe and CISe. We'll see that the such behaviour will be recovered when passing to the mBJ XC potential.

The "optical" band gap between the occupied (blue) and vacant (orange) states is negligible (98 meV for CGSe and 2 meV for CISe) yet existent in Γ. The experimental values (to be discussed below) are much larger.

The contribution of different orbital states in different bands and k-points is shown in Fig. 3.5, for the case of CISe. The figures of similar design will also come about in the following, so some explanations can be made here. The solution of the Kohn-Sham equation yield eigenvalues, i.e., the band energies E(k), and eigenvectors. The eigenvector reflects the participating of different basis functions at the given solution. The basis functions of the WIEN2k code are augmented plane waves, i.e., the (quite numerous) plane waves passing (within the atomic spheres) into the numerical functions

Band structures of pristine systems calculated with GGA and mBJ

centered at atoms and numbered by the atomic site and the (l, m) (orbital, magnetic) quantum numbers. Within a given eigenvector, one can sum up the squares of all its components which refer to the given atom and the given l index, whatever the m values and the plane waves participating. The corresponding sum is a measure of the given (atom, l-value) participants in the eigenvector concerned, and is plotted as a circle of proportional radius around the corresponding eigenvalue. The general scaling factor (for the whole plot) is chosen out of graphic convenience in each plot independently. The underlying band structure (neglecting the circle radii) is the same in all three panels of Fig. 3.5 (and also the same as in the right panel of Fig. 3.4). As was argued before, the bands at -13 • • • -14 eV are of almost exclusively Se4s character (as is evident from the middle panel), and otherwise a noticeable Se s-contribution is only seen in the band #60 (the third one back from the top of the valence band) near Γ, where this band ascends to -0.4 eV. In the left panel one sees that this band has also a noticeable In s contribution near Γ. Otherwise the In s participation is shared between the occupied bands at -7 • • • -6 eV and the vacant bands at 0 • • • 2 eV. From the right panel one sees that the Se p states contribute into a large number of bands, from -7 eV upwards (including the "nominally In4s" band at -7• • • -6 eV).

mBJ

The use of mBJ XC potential (explained in Subsec. 1.2.3) modifies the band structure, the most important effect being the general shift of the vacant states upwards, that results in an increase of the band gap (or, sometimes in opening of the gap in the systems which failed to be correctly identified in semiconductors or insulators in a calculation done within LDA or GGA). Fig. 3.6 shows this modification, to be compared with the GGA band structure in Fig. 3.5. What happens is a slight general shift of the unoccupied bands upward, in addition to which somehow stronger perturbations happen at both sides (within about ±5 eV) from the "Fermi energy" (separating occupied and vacant states).

The band gap opens to 0.46 eV in CISe. Interestingly, the situation with band dispersions in the vicinity of the optical gap at Γ is now "healed" and becomes similar to that in CGSe: the upper occupied bands, #60 to #62, ascend towards Γ and become degenerate in it, whereas the In s character is completely pushed from these bands into the vacant ones, #63 and #64. Comparing the left panels of Figs. 3.5 and 3.6, one can get an idea of what is going on: in the GGA calculation the lower vacant band, with pronounced In s character (and also with an important Se p contribution) would in fact dive by ∼ 0.4 eV into the valence band. However, its crossing with the convex upper valence band #62 is avoided, and the band structure arranges itself with a "technically zero" band gap, pushing the occupied band #60 downwards in Γ. A passage to mBJ shifts the conduction band generally upwards by ∼1 eV; the "immersed" lower vacant band #63 then "gets free" from the valence band and stabilizes in the "qualitatively correct" configuration of bands with its minimum at 0.46 eV. It is remarkable that the effective upward shift of the lower conduction band in CGSe is of about the same net value as in CISe. The difference is that in CGSe a tiny band gap (0.08 eV) existed already with the GGA, and was increased with mBJ up to 1.06 eV. The difference between the band gap values in CGSe and CISe is thus correctly established, even as a systematic underestimation of experiment data by about 0.5 eV still takes place. The discussion about the band gap and its variation with concentration follows in the next subsection.

The comparison of band structures of CGSe and CISe as calculated with mBJ is shown in Fig. 3.7; this figure has to be compared with Fig. 3.4 to see the effect of mBJ vs GGA-PBEsol. Besides the opening of the band gaps, a slight contraction of the width of the upper "blue" bands (with large Se4p and Cu3d contribution) can be noted. This effect is also seen in the increase of the density of states (scaled identically with that in the DOS panel of Fig. 3.4).

For reference purposes, the density of states of CuInSe 2 is shown in Fig. 3.8; that for CuGaSe 2 would be very similar. Note that the partial DOS are defined from inte- and partial (per atom, the lower panels) densities of states in CISe, as calculated by WIEN2k with the mBJ XC potential.

Mixed systems; band gaps varying with concentration

states give rise to two separated bands situated at ∼ -7 • • • -6 eV. This reflects the fact that the Ga4s state is slightly deeper than In5s. In larger supercells Cu 4 In 3 GaSe 8 and Cu 4 InGa 3 Se 8 , where the total number of bands doubles, the structure of these in this region changes correspondingly: in the first case, a deeper single band is separated from three placed ones above it whereas in the second case, the isolated band is higher. The Ga3d / In4d states, of which Ga3d is slightly deeper, are dominating in the narrow band at ∼ -16 eV, (in fact, a bundle of bunds, five from every In/Ga atom participating) which acquires a structure according to In/Ga concentration: its weight its shifted the upper edge in CuIn 0.75 Ga 0.25 Se 2 , to the lower edge in CuIn 0.75 Ga 0.27 Se 2 and evenly distributed in CuIn 0.5 Ga 0.7 Se 2 . As was discussed in Chapter 2, the 1:3 and 3:1 supercells imply the doubling of the unit cell (from body-centered tetragonal to simple tetragonal) and shrinking of the BZ. The Z-Γ path of the bct BZ, when mapped onto the single tetragonal BZ, starts at Γ of the adjacent BZ and hits the BZ surface at half-way. The corresponding labelling is used in the k-path of the first and the third panels of Fig. 3.9, even if "physically" the path is the same in all three plots. The doubling of the unit cell leads to doubling of the number of bands, which are "reflected" at the zone boundary at "new" Z and folded back into the BZ, like in the textbook example of doubling the cell in a linear chain.

It is remarkable that the vicinity of the band gap (at Γ) is only slightly quantitatively affected, compared to cases of pure compounds. In all cases, independently on the total number of bands and their folding over the BZ, the band structure in the immediate Chapter 3 vicinity of the band gap is cleanly made by a single (the lowest) among the vacant bands, protruding downwards at Γ, and three highest bands among the valence ones, which protrude upwards in Γ and are degenerate in this point. This consistency justifies somehow the choice of ordered supercells for simulation of solid solutions. One notes a regular increase of the band gap with Ga concentration. This can be seen in Fig. 3.10, in comparion with experimental results extracted from [START_REF] Schnohr | Compound semiconductor alloys: From atomic-scale structure to bandgap bowing[END_REF]. The magnitudes are systematically underestimated by ∼ 0.5 eV, whereas the slope (the total variation of the band gap from CISe to CGSe) of ∼ 0.6 eV is very close to experiment. One can speculate even that the band gap bowing of b=0.2 eV, according to [START_REF] Schnohr | Compound semiconductor alloys: From atomic-scale structure to bandgap bowing[END_REF], is grasped correctly; the upward kick at x=0.5 may well be an artefact of having a too ordered supercell; any disorder would tend to reduce the band gap value.

Summary on chalcopyrite-type structures

• According to WIEN2k calculations, PBEsol-GGA provides a very accurate description of lattice parameters in case of CuInSe 2 , but slightly (by ∼0.1%) overestimates the c parameter in CuGaSe 2 .

• The band gap is direct at all concentrations; it occurs at Γ between the triply degenerate highest occupied bands and the isolated lowest vacant band.

• The band gap values are systematically underestimated by ∼ 0.5 eV, but the slope (and, probably, bowing) in dependence on concentration is faithfully reproduced. It seems that the GGA calculation for CGSe would show in fact a "physical" overlap of the highest occupied and the lowest unoccupied bands. Due to an avoided crossing of bands, the system "technically" opens a tiny band gap and remains semiconductor.

Chapter 4

Results for GaSe and InSe systems

Introduction

The main aims of the study of the GaSe, InSe crystals, like in relation to previously addressed chalcopyrite-based systems, were the accurate lattice relaxation (attributing the major attention to the performance of PBEsol XC functional), and the description of the band structures and the band gap (especially the performance of the mBJ XC potential). The peculiar features of the compounds under study are : (i) the layer structure in which the double anion-cation-cation-anion layers are separated by "van der Waals gap" and, presumably, an importance of dispersion interaction for fixing the interlayer distance, hence problems for "conventional" DFT methods;

(ii) the presence of competing polytypes with, presumably, very close energies. For (i), the VASP calculations have been organised in parallel with the WIEN2k ones, because the available version of the VASP code takes into account (in parametrized form) the dispersion interactions within the Grimme formalism (see Chapter 1). For (ii), special efforts have been applied to arrive at utmost precision in order to resolve the energy / volume curves for different polytypes. Like in the previous chapter, the outline of lattice relaxations and related observations (for both GaSe and InSe systems) will be followed by the discussion on the structures, making reference to the copper-containing chalcopyrite-type structures and explaining the effect of mBJ.

Optimized crystal structures

As was mentioned in Chapter 1, the structure relaxation in WIEN2k and VASP codes was performed in technically different way. In VASP, the procedure is more automatized, whereby one can choose between different relaxation regimes, freezing or releasing either lattice parameters or internal coordinates. The WIEN2k way is more rigid and more Chapter 4 dependent on human intervention; the sequence of fittings (and minimisations) of the total energy as function of volume and c/a is alternated by "sessions" of optimizing the internal coordinates (by minimizing forces on all atoms) at fixed lattice parameters. The results of this analysis for different polytypes, after both WIEN2k and VASP and in comparison with experimental data, are shown in the following figures and tables.

Lattice parameters; an overview of different XC schemes similarly small (minimal) lattice parameter in the hexagonal plane, a, but different number of double layers. For the sake of comparison, the c values are "normalized" in the following to the number of structural units per unit cell (two for β and ε, three for γ, four for δ. This will be further on applied also to the discussion of total energies). Moreover it can be mentioned that the WIEN2k calculations for the γ phase have been technically done for the rhombohedral unit cell whereas the VASP calculations used the hexagonal unit cell, tripled with respect to the rhombohedral one. The experimental results are only available for several phases (β,ε-GaSe and β,γ-InSe) studied at different conditions; however they reveal the closeness of a, c parameters throughout the compound given (as well as the closeness of results from different experiments). The scattering of the calculation results, also those obtained within the present work with different XC schemes, occurs at much larger scale. First it should be noted that the GGA-PBE results ("original" PBE, see Chapter 1) are very far off, considerably overestimating both a and c parameters. This can be "repaired" by switching to the PBEsol scheme (applied within the WIEN2k code), or by including dispersive interaction on top of PBE within one or another variation of the Grimme scheme (indicated PBE+D2, PBE+D3 or PBE+D3-BJ) implemented in the VASP code. (Most of the trial calculations with VASP have been done by Michaël Badawi to whom I owe the completeness of these tests). It is obvious that a combination of PBEsol with the Grimme correction (PBEsol+D3, PBEsol-D3-BJ), also tried with VASP, too much overbinds in all cases. As a result of practical observation rather than of a profound analysis, the PBEsol (with WIEN2k) and PBE+D2 (with VASP) are retained for the following detailed inspection of band structures, enery/volume curves, etc.

A more attentive look at the results of Fig. 4.1 reveals that, whereas the situation with GaSe can be considered as satisfactory (good agreement of PBEsol and PBE-D2 results between themselves and with the experiments available), for InSe it is by far not the case. PBE+D2, while correctly predicting c value, yields an apparently too short (by ∼1%) in-plane lattice constant; PBEsol is performing much better in this sense. The PBEsol, in its turn, has a strange problem with -of all phases -ε-InSe, for which Chapter 4 it predicts the c parameter longer as expected by as much as 6%. This seems unlikely to have a genuine physical reason; however, a careful search for a technical error in repeated calculations, attentive comparisons with those for the β phase etc. so far did not produce any evidence of error. In the lack of explanation, the results for the ε-InSe phase have been withhold in the publications but are exposed here as they are.

A vague interpretation of the difference between GaSe and InSe, in what regards the use of PBE+D2, could be a presumed larger "degree of covalence" and delocalization of

Optimized crystal structures

electron states in the latter. So to say, the double layers in GaSe are more "isolated" by the vdW gap and hence "underbound", so an explicit inclusion of dispersion interactions via the Grimme's corrections is needed to fix the situation. In InSe, the situation is already "covalent enough", so that an inclusion of D2 on top of PBE (or, even worse, on top of PBEsol) overbinds too much. It is remarkable that PBEsol does a fair job in all cases (at the exception of "mysterious" ε-InSe), thus providing a reasonable practical compromise.

The results of thus "selected" calculations, i.e., WIEN2k with PBEsol and VASP with PBE+D2, are further specified in Table 4.1. According to them, the net thickness of the double layer h (i.e., Se-Se distance across the layer) and the vdW gap d somehow vary over the polytypes (as a manifestation of different stacking), whereas the (intralayer) bond lengths l remain quite stable. Referring to the above mentioned anomaly in the optimized c parameter of the InSe-ε phase, one can note that it is namely the interlayer distance (the vdW gap) that is largely overestimated (∼3.5 Å vs ∼3.0 Å for the other phases), not the distances within the layer.

It can be emphasized again that the PBEsol calculations yield very accurate estimates (vs. experiments) for the a lattice parameter, whereas the c/n is generally Table 4.2: Interatomic distances in GaSe and InSe polytypes according to WIEN2k PBEsol calculations, separately for cation-anion nearest-neighbours and next-nearest neighbours, cationcation bonds within the double layer (perpendicular to it), anion-anion interlayer (connecting Se atoms across the vdW gap), anion-anion [001] across the double layer (the layer thickness) and anion-anion in-plane (the a lattice constant). For the δ phases, all the available entries in each category are shown. overestimated by ∼0.1 Å. This overestimation is quite accurately corrected for β, γ and δ of GaSe via inclusion of the Grimme correction (PBE+D2). A more careful insight reveals that the vdW gap d comes out a bit too short, that is compensated by a bit overestimated double layer thickness h. This "deformation" hardly affect the Ga-Se or In-Se bond lengths l which are almost identical after PBEsol and PBE+D2. We will see some similarities in their band structures as well.

cat-ani

To conclude the section about equilibrium structures, the resulting interatomic distances are shown in Table 4.2. They depend on the internal coordinates of atoms but are often explicitly given in experimental or theory papers. For the δ phase, there are several nonequivalent positions and hence different "species" of various interatomic distances, all of which are listed in the Table 4.2. Even as the scattering of these values is very small, one can note, à toutes fins utiles, that the largest disparity is in the anion-anion interlayer distances, hence in the width of the vdW gaps, and not in the thickness of double layers.

Band structures and Densities of States

The band structures discussed calculated under different conditions and shown in this section will share much similarity. However, we will try to concentrate on their small yet important differences. The general composition of band structure, from individual electronic states of participating atoms, is already known from the previous chapter, only that now there are no Cu3d states participating. The other difference is, the compounds in question do not satisfy matching the formal charge, e.g., Cu + Ga 3+ Se -2 2 , so we deal with a partially occupied Se4p shell, and the band gap (as there will be a band gap!) develops between bands with Se4p participation. Apart from this, the semicore Ga3d / In4d states, as well as Se4s and the "below the valence band" Ga4s / In5s narrow yet dispersing bands can be found at their expected locations.

An analysis of orbital compositions in several symmetry points of the BZ has been done already by [START_REF] Robertson | Electronic structure of GaSe, GaS, InSe and GaTe[END_REF], in relation with the latter's calculation within a tight-binding model. Later on, [START_REF] Gauthier | Band gap variation of gallium selenide under high pressure[END_REF] suggested a simple scheme that is reproduced in Fig. 4.2, which has been further elaborated in the work by Gauthier et al. [1989] -see Fig. 11, 12 therein. Gauthier et al. argued that the splitting at Γ marked E (1) in the figure is due to Ga-Ga (4s) intralayer interaction, the E ′ (1) splitting in M -due to Ga 4p xy intralayer interactions, whereas the "secondary" splittings E BV and E BC are induced by inter layer interactions involving mostly Se 4p z states, and E ′ (1) BC is due to interaction in M between Ga 4p z states. These simple considerations generally remain valid in view of later calculations, including the figures presented below in this chapter.

A subtle issue is the band foldings, on passing from smaller-cell to larger-cell poly-Th- -E 1 splitting of the s-p gallium le v els at the r point.

-rOE .J:• E' "I !Ed -[ .-•-1 1 rc-•-_t __ : ~ ! 1 1 r M

Band structures and Densities of States

-E 1 < 1 ) splitting of the P KY gallium levels at the M point.

-Env• EBC splLtting of the P selenium levels in the vilence and conduction bands.

-E' sc splitting of the P 2 gallium tevels in the conduction band. types. The related splittings are the most visible along the (short) r-A k-path in the [001] direction. Finally, the use of the mBJ XC potential shifts the bands and affects the electronic states more in one parts of the BZ than in the others. We start with an outline of a "naiven model of band dispersion in a hexagonal double layer. As a.

four sites / four states model, it neglects the "richness" of the hybridization between the Gap and Sep states, but has an advantage of being exactly solvable. Further on, a number of band structures is shown. The comparison will be concentrated on the cases [the same XC potential, different polytypesJ, (a given phase, WIEN2k vs VASP] and [a given phase, mBJ vs PBEsol]. The +D2, +D3 corrections aimed at including the dispersion interactions do affect only the total energies / forces and have no immediate effect onto the band structures ( only via the modified relaxed geometry), so they are excluded from this consideration.

Band structures with PBEsol calculated by WIEN2k

Band dispersions for different polytypes, as calculated with P BEsol, are shown in Fig. 4.3. They can be understood as t hose of a ((prototype" double layer, upfolded according to the number of units (double layers) per unit cell and slightly modified in the course of structure relaxation in each particular system. As is well known from the textbook example of dispersion of electronic states on a linear chain, and the effect of doubling the unit of such chain, a given state in r is split in two , at the corresponding ascending and descending dispersion branches get degenerate at the zone boundary. The doubling (or tripling, quadrupling) the cell occurs here along the [001] direction, and the zone boundary is at A. Accordingly, the band structure of ( doubled) f3 and e PBEsol XC functional and the conventional PBE. In view of these tiny differences, the agreement looks astonishingly good, especially taking into account that the "inner mechanics" of the two calculation methods are quite different. For one thing, the Se3d states counted within the valence band in WIEN2k are omitted from the valence band (and managed within the projected augmented wave technique) with VASP. One can see that some band splittings at Γ are just a bit larger in the VASP calculation than with WIEN2k. Being generally satisfied with the agreement, we did not further explored whether these tiny differences can be traced to the crystal structure, or XC functional, or the calculation method (basis etc.) used.

Band structures and Densities of States

Band structures and Densities of States

The band structures of InSe phases are quite similar; an explicit comparison of GaSe to InSe will be done for the case of calculations performed with the mBJ XC potential.

Band structures calculated with mBJ; band gaps

The set of band structure plots for the GaSe polytypes, as calculated with mBJ, is shown in Fig. 4.6, and their InSe counterparts -in Fig. 4.7. As could have been The numerical values of band gaps from the present calculations and experiments are collected in Table 4.3. Whereas, as expected, the band gap is increased while making a calculation with mBJ, some more specific observations can be formulated: (i) All the phases are semiconductors with appreciable band gaps already according to "conventional GGA". The "augmentation" of the GGA gap due to mBJ (by ∼1 eV) is much larger than in CGSe, CISe (∼0.6 eV). (ii) GaSe reveals a much larger (by ∼0.8 eV) band gap than InSe, that is consistent with the situation in chalcopyrite-type systems (the gap in CGSe in by ∼0.6 eV larger than in CISe); this trend exists with GGA and is "translated" into the mBJ results. (iii) According to PBEsol, GaSe has a direct Γ-Γ gap, that is replaced by a indirect one, since the conduction band acquires a minimum in M slightly deeper than in Γ. In InSe, no such "transformation" occurs, the band gap remains direct in all cases.

Band structures and Densities of States

(iv) The nature of the band gap (indirect in GaSe, direct in InSe) is correctly predicted by calculation, and the numerical values are in unusually good agreement with experiment. Chapter 4

Scans of total energy vs volume; the relative stability of phases

The calculated equilibrium structure gives no clue as for how stable, or unstable, the particular phase is compared to possible competing ones. In fact, the static energies of these competing phases allow to judge about their relative stability only in a very indirect way. A more systematic study consists in inspecting the variation of the total energy depending on some structure parameter, typically -the volume (or, pressure). This helps to figure out how steep (or flat) the local energy minimum is and how easily the system can be, under certain conditions, driven into another crystal structure. Such tests have been done for all phases studied of both GaSe and InSe systems, using both WIEN2k and VASP methods. The idea of using the both method was to get the general confidence in the accuracy in a procedure where a quite high precision was expected to be needed in order to resolve the presumably very similar behaviour of different polytypes. The other reason was to make use of the option, provided by the VASP code, to perform full relaxation of the crystallographic parameters (c/a and the internal coordinates) for every trial value of volume. This is, in principle, a correct method to study the Energy(Volume) or Energy(Pressure) relation, that would allow to extract the bulk modulus (if needed), or other elastic parameters, applying corresponding constraints. The WIEN2k code, from its side, is less flexible in performing volume-dependent scans, because anything other than the simplest uniform expansion/compression (controlled by the single lattice parameter as the universal unit length), becomes technically quite complicated.

As an attempt of a meaningful comparison between the results by VASP and WIEN2k, Fig. 4.9 depicts the scans of total energy as function of unit cell volume, which pass through the fully optimized structure for each respective polytype and explore the uniform expansion / compression of the unit cell from this optimized geometry. The absolute total energy values do not have any special meaning, but the relative energy scale make sense and should be comparable between the methods. Indeed, the general shape of the nearly parabolic scan in the right panel depicting the WIEN2k results quite perfectly coincides with the steep parabola in the left panel, marked "c/a fixed", which reveals the same type of distortion applied in VASP calculation. Unfortunately, the agreement between the methods in what regards fine details, i.e., the placement of individual parabolas, is not that perfect. The less trivial observations can be systematized as follows:

-According to both calculations, the β phase is the ground-state one, producing the lowest-lying parabola.

-The spread in energy between the polytypes, for a given volume, is very closearound 1 meV (per two formula units, f.u.) in VASP and around 1.5 meV according to WIEN2k. for different phases of GaSe, using the PBEsol XC potential. The steeper set of curves (those for WIEN2k, and those labelled "c/a fixed" for VASP) correspond to the uniform scaling of the structure at equilibrum. The softer set of curves (labelled "c/a relaxed" for VASP) correspond to the fully relaxed structure for every volume.

Total energy vs volume for polytypes of GaSe and InSe

-The equilibrium volume, not much different throughout the polytypes, is ∼100 Å3 (per two f.u.) according to VASP and ∼98.5 -99 Å3 according to WIEN2k.

-The placement of the curve corresponding to the γ polytype shows the largest difference between the two calculations: it is the highest one according to VASP but almost competing to the ground-state β curve, according to WIEN2k. In this relation it should be mentioned that the γ phase was treated as rhombohedral one in WIEN2k but as a hexagonal one (with triple amount of atoms in the unit cell) in VASP, so that the systematic error is probably the largest in the case of this phase.

-The ε and the δ phases are consistently identified as the highest one and the intermediate one, correspondingly, according to both calculations. Taken together, this seems to reveal a comforting level of agreement, considering a very different background of calculation methods and the delicacy of the differences to identify. We tend to look at it as a justification of VASP results by comparison with a, a priori, superior method in what regards the precise total energy calculations for a given geometry, whereas the relaxed structure data might well come out more accurate in VASP, due to more sophisticated and less "human-dependent" algorithm implemented earlier, [START_REF] Peressi | Structural and electronic properties of Ga 2 Se 3[END_REF] studied, in a first-principles calculation, an orthorhombic body-centered structure and predicted this material to be semiconducting, with a direct gap of 1.26 eV already in LDA (i.e., probably strongly underestimated).

Relaxed crystal structures from calculation

To maintain a link with the previously discussed results on 1:1 systems, in this chapter only hexagonal structure model will be discussed, which might seem somehow hypothetical for Ga 2 Se 3 . However, we note that the monoclinic phase is only moderately distorted, as compared to hexagonal, and was even mistaken for the latter is some early works. Moreover, the study of the driving force of the monoclinic distortion in Ga 2 Se 2 , as compared to the lack of such tendency in otherwise apparently similar In 2 Se 3 , might be enlightening for future work. Another interesting issue is that, independently of the monoclinic distortion, the underlying structure of Ga 2 Se 3 seem to be definitely "zincblende-like", whereas different phases of In 2 Se 3 fall into "wurtzite-like" category. The analysis of this discrepancy, might also deserve attention in subsequent studies.

For plausible structure models used in the present calculations, we follow the work by Ye et al. [1998], which explains different ways of ordering the vacancies in the underlying wurtzite structure, the "vacancies ordered in screw form" (VOSF) and the "layer structure phase" (LSP). These models have been explained with some details in Chapter 2.

The VOSF unit cell is a √ 3× √ 3×3 multiplication of the wurtzite one and contains 9×4-6 = 30 atoms, hence 6 formula units. The corresponding space group can be either P 6 1 (Nr. 169) or P 6 5 (Nr. 170), according to the screw being left-hand or right-hand. The LSP possess the wurtzite lattice vectors in the basal plane (without multiplication), and the vacancies are assembled into missing cation layers, so that instead of "double layers" (VI)-(III)-(III)-(VI) the system consists of "triple layer" packages (VI)-(III)-(VI)-(III)-(VI), separated by van der Waals gaps between the opposing anion layers. The structure, with its possible variations, is explained by Fig. 10 of [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF].

Whereas such triple layers can be assembled in an infinite number of possible sequences, as is generally the case with hexagonal polytypes, [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF] distinguishes between the two simplest stackings, according to whether the Se sheets within each triple layer follow the "wurtzite-like" sequence, A-B-A, or "zincblende-like" one, A-B-C. Taken together, this yields two possibilities, shown in Fig. 10

Relaxed crystal structures from calculation

The structure of VOSF phase is shown in the left panel of Fig. 5.1, to explain the positions of 2-coordinated and 3-coordinated Se atoms, to be referred to later. The right panel of Fig. 5.1 depicts the layered phase.

For all three structure models, the lattice parameters have been optimized in a sequence of WIEN2k calculation steps (this involved also optimization of internal coordi-

Band dispersions and orbital contributions

DFT calculation (also in those done with PBEsol), seem to be stronger underestimated in Ga-Se systems than in In-Se ones. Consequently, the Ga 2 Se 3 lattice looses a bit on insertion of "vacancy layers", whereas the In 2 Se 3 lattice becomes a bit tighter due to a stronger covalent part in Se-Se interaction across the vdW gap.

The study of the relative stability of polytypes would require to scan there energy/volume curves near the corresponding equilibria, as was done in Chapter 4 for the 1:1 phases. This was not concluded at this moment.

Band dispersions and orbital contributions 5.3.1 Layered phases

Band dispersions of two layered phases and of the VOSF, as calculated within GGA (using PBEsol XC functional), are shown in Fig. 5.2. All these phases are hexagonal; the labelling of high-symmetry points in the BZ can be found in Fig. 2.6. However, the dimensions of the BZ and the number of bands are different: VOSF has roughly √ 3 times longer in-plane lattice parameter and 2/3 times the c parameter than the LSP; Note that the layered phases are not semiconducting according to GGA.

Chapter 5 correspondingly, the number of bands is twice that of LSP. The band structures in Fig. 5.2 are aligned by the "Fermi level" (chosen as zero at the energy scale), that is, for the VOSF phase which exhibits a band gap, by the top of the valence band. This has an inconvenience that well localized In4d states, as well as the Se4s band, are not aligned. Concerning the Se4s states, one notes a narrow non dispersive band (at ∼ -11.5 eV in 5.2) which is quite separated (by ∼2 eV) from the "main" (dispersive) Se4s band (at -12 --14 eV). This split-off state is due to the 4s states of Se atoms which are in incomplete coordination and therefore miss to be efficiently hybridized with the states of neighbouring atoms. In the layered structures, these missing-coordination atoms are the Se atoms at the surface of the Se-In-Se-In-Se layer, or, more specifically (since the layer is not symmetric), the "lonely" Se atoms which are on top of the cation beneath, and bonded only with the later.

To remind the geometry of the layer, we refer to the right panel of Fig. 5.1 which is labelled Ga 2 Se 3 ; the structure of In 2 Se 3 would be identical. The "lonely" atoms in question are those labelled as "Se (top)" in the figure cited. The Se atoms on the

Band dispersions and orbital contributions

opposite surface of the layer (marked "bottom") also have an incomplete coordination, but of much less dramatic type: they bond to three cations in the same (warped) wurtzite-like plane; only the fourth their bond is not saturated. The Se atoms in the middle of the layer have full tetrahedral environment. This different environment and different saturation of covalent bonds for different Se atoms manifests itself also in the density of states of correspondng atoms, shown in Fig. 5.3 (for the case of Ga 2 Se 3 ; the trends are very similar for In 2 Se 3 -LSP). The Se (top) atoms not only have very narrow and less deep 4s state, but also their Se4p band is relatively narrow and shifted upwards on the energy scale, compared to other Se atoms. In principle, the upward shift of 4s and 4p can be anticipated, as the crystalline potential at the peripheric "top" atom is less deep than within the layer.

Coming back to Fig. 5.2 and the absence of band gap for the LSP phases, one can see that this is due to an overlap, that happens only in some vicinity of the BZ center (including Γ and A) of convex upper valence bands and protruding down lower conduction bands, similarly to what we observed in CuInSe 2 according to GGA calculation.

Vacancy-ordered phase

The band structure of the VOSF phase (Fig. 5.2) also contains a narrow split-off band at ∼ -11.5 eV, and its origin can also be traced to a reduced coordination of a Se atom. Indeed, the crystal structure of VOSF, as was explained in Chapter 2 and is shown in the left panel of Fig. 5.1, includes 3-coordinated and 2-coordinated Se atoms. Their local densities of states (Fig. 5.4) and contributions to the eigenvectors in different parts of the energy dispersion pattern (Fig. 5.5) can be analyzed separately. From Fig. 5.5 that shows orbital contributions from 4s and 4p states of three-coordinated and twocoordinated Se atoms it becomes evident that the split-off band among those related to the Se4s states is brought about almost exclusively by the states of the 2-coordinated Se atoms. Correspondingly, the Se4p states of these atoms, within the valence band, are dominating towards the top of the valence band. These trends can also be seen in plot of the densities of states, shown in Fig. 5.4.

Band gaps

The GGA calculation predicts the band gap of 0.916 eV for Ga 2 Se 3 and 0.552 eV for In 2 Se 3 , both being in the VOSF phase. Such relation is consistent with qualitative trends found in other earlier discussed In-based and Ga-based compounds. Layered phases, according to GGA, have no band gap. With mBJ XC potential, the band gap predictably increases to ∼1.986 eV for VOSF-Ga 2 Se 3 and ∼1.587 eV for VOSF-In 2 Se 3 . In 2 Se 3 (PBEsol) 3-coordinated Se 2-coordinated Se Γ (in Ga 2 Se 3 , not so markedly in In 2 Se 3 ).

Band dispersions and orbital contributions

The experimental works report a number of estimates of the band gap of In 2 Se 3 , which are sometimes difficult to bring in relation with the crystal structures, due to ambiguity in the latters' attribution (as seen in Table 5.1). [START_REF] Julien | Optical properties of In 2 Se 3 phases[END_REF] cite the optical gap value of 1.812 eV (at room temperature) or 2.05 eV (extrapolated to zero temperature) for the phase they refer to as γ-In 2 O 3 , a high-temperature one (cf. "van Landuyt γ" in Tab. 5.1) but which presumably can be stabilized also at room temperature. Further on, [START_REF] Julien | Optical properties of In 2 Se 3 phases[END_REF] cite the optical gap of the α phase as ∼1.35 eV at room temperature, to become 1.560 eV when extrapolated to zero temperature. In principle, both "Likforman α" and "van Landuyt γ" from Table 5.1 could be reasonably associated with the VOSF, based on the similarity of their respective lattice parameters. Under this angle, one can conclude that the mBJ calculation yield the band gap either (assuming α=VOSF) in unexpectably good agreement with exper- iment, or (assuming γ=VOSF) slightly underestimated, as is typically the case in many other semiconducting systems.

For Ga 2 Se 3 , the situation is more delicate, because the most relevant phases under discussion seem to be not of hexagonal structure. For instance, [START_REF] Takumi | X-ray, Raman and photoluminescence study of vacancy ordered β-Ga 2 Se 3 under high pressure[END_REF] report the fundamental absorption edge of 2.3 eV (at ambient pressure, among their studies done at elevated pressures) for the (monoclinic) β-Ga 2 Se 3 phase, whose structure was specified by [START_REF] Lübbers | The crystal structure of β-Ga 2 Se 3[END_REF]. As is seen from Fig. 2 of [START_REF] Lübbers | The crystal structure of β-Ga 2 Se 3[END_REF], the nearest-neighbour arrangement and the distribution of cation vacancies have similarities to the correspondent characteristics for the VOSF, only that the underlying structure of Ca 2 Se 3 is zincblende and not wurtzite. Similarly, [START_REF] Okamoto | Anomalous anisotropy in the absorption coefficient of vacancy-ordered Ga 2 Se 3[END_REF] refer to zincblende-related phase of Ga 2 Se 3 in which Ga vacancies occupy every third (011) or (01 1) plane, and estimate the band gap to be 2.1 eV. Such reasonable agreement with the band gap value calculated for VOSF Ga 2 Se 3 (as well as the observation that the band gap exceeds that of In 2 Se 3 by ∼0.4 eV) may not be fully fortuitous.

Calculations with mBJ for layered phases encountered some technical problems (very unstable start-off while initializing the mBJ procedure). Comparing the "extrapolated" value of increasing the gap due to mBJ, say ∼1 eV, with the "negative gap" in Fig. 5.2 (the overlap of the valence and conduction bands), which is of about the same order, 5.4. Assessment of results on (III) 2 Se 3 systems one can conclude that the opening or not of the gap is possible but not granted.

Assessment of results on (III) 2 Se 3 systems

The work done so far on In 2 Se 3 and Ga 2 Se 3 leads to the following conclusions.

1. √ 3× √ 3×3 superstructures on the basis of wurtzite have been very actively under discussion in relation with In 2 Se 3 . The VOSF phase is a plausible realisation of such model. The ordered distribution of cation vacancies imposes the presence of 2-coordinated and 3-coordinated Se atoms, which quite markedly differ by their bond lengths to neighbouring In atoms, and give rise to a very pronounced features in the band structure / densities of states (narrow peaks hinting for localized states). According to DFT calculations, this phase is semiconducting already when using the GGA. The use of the mBJ XC potential opens the band gap even further. In fact the predicted band gap values obtained with mBJ are in quite good quantitative agreement with experiment, but this should be looked at with caution, for two reasons. First, the reference In 2 Se 3 phases on which the experiments have been done may only presumably be related to VOSF. Second, the Ga 2 Se 3 phases used in experiments are zincblendebased and may, at best, only provide the nearest coordination and distribution of cation vacancies similar to those in VOSF phase, but not the long-range order nor the local symmetry of the later.

2. The layer structures, whose in-plane lattice constant is essentially that of wurtzite, make the other "family" of phases involved in many discussions. Obviously, the arrangement of cation/anion layers, within the given stoichiometry, can be realized by many possible ways. The LSP1 and LSP2 are just two models which seem plausible, in view of existing literature. First-principles calculations for them encounter a problem that the band structure of these materials, according to GGA, is not semiconducting. By itself, it does not necessarily "devaluate" the predictions done in the calculations concerning the equilibrium structures. However, the necessary caution should be applied in assessing these data. It is likely that applying the mBJ treatment will result in opening a band gap. Unfortunately, this won't be necessary useful for "improving" the structure predictions, because the mBJ XC potential was designed just in view to "correct" the band energies and does not include any "mechanism" to correct at the same time the total energies and forces acting on atoms.

3. In any case, on the base of the results so far obtained, one can conclude that the reduced coordination of some Se atoms in the layer structures manifests itself in somehow similar way as in VOSF. Namely, the lengths of bond connecting such lowcoordinated anions to In atoms are very noticeably reduced. Again in similarity with what happened in VOSF phases, localized electronic states have been formed related to these low-coordinated atoms, characterized by narrow peaks in the density of states. Chapter 5

4.

In what regards Ga 2 Se 3 , the strictly hexagonal models like VOSF or LSP seem to be too idealized, because numerous experimental indications exists towards distortion into monoclinic structure. However, in the present work the idea was to trace the Ga-based and In-based selenide systems in comparison, through different compositions and phases. An analysis of instability versus changing from hexagonal into monoclinic structure may become an interesting option for future first-principles calculations.

Conclusion

This works brings some studies to supposed termination, but quite a number of issues failed to find their due resolution. The simulations on chalcopyrite-type structures were confined to somehow too ordered model systems. Various III-VI polytypes are predicted to have astonishingly close energies, but only few of these systems are realized in nature; the reason for this remains obscure. The effects of lattice vibrations and of spin-orbit interactions were completely left out of the scope of this work, yet they might play important role in some situations. The colleagues doing experiments in vibrational or electron spectroscopies have big expectations from theory but they can not always be rewarded... This work signifies a long road which I passed since first confronting with the issues of first-principles calculations for materials, through many small discoveries, victories, annoyances and despairs towards a stage where some parts of work are concluded, the others ought to come to conclusion very soon, and yet others did not went as expected, and demanded to take new turns and to change the plan on the road. I learned about semiconductors and photovoltaics, crystallography and phase diagrams, calculation methods and their pitfalls. It was rewarding to learn how first-principles calculations are able to do very accurate predictions, and what various information can be extracted from very basic electronic properties. There are yet so many interesting problems to look into, and so many interesting things to learn... This work resulted from my PhD years offered me various and many defiances regarding the study of semiconductors of binary, ternary, big crytals and alloys. It offered me many intellectual and practical skills leading me to a new way of thinking scientifically and logically... And a new proficiency and expertise to solve the scientific problems more effectively.

It was not originally intended, but utimately turned good, to apply two methods of first-principles calculations, both considered as very reliable, in comparison. Done out of frustration about particular technical problems, it turned out very useful to give some confidence in particularly sensitive calculations, and to estimate the real accuracy and credibility offered by state-of-art theory.

Figure 1 . 1 :

 11 Figure 1.1: A fragment of the .struct file for WIEN2k calculation of the ε phase of GaSe.

Figure 1 .

 1 Figure 1.3: A k-convergence test for different phases of GaSe, after WIEN2k calculations.

Figure 2 . 1 :

 21 Figure 2.1: Tetragonal unit cell of CuInSe 2 in the chalcopyrite structure in a perspective view (left) and in the top view (right).

2. 1 .

 1 General introductionTable2.1: Wyckoff positions in chalcopyrite-type CISe or CGSe (space group I 42d). The missing coordinates in the tetragonal cell (Fig.2.1) can be obtained by the translation (

4 (

 4 Nr 82), with Cu atoms in (a) and (c) positions, In in (b) / Ga in (d), or other way around, and Se in the (g) positions. Replacing one of (so far equivalent) Ga atoms by In would result in In:Ga=3:1 supercell (shown in the right panel of Fig. 2.2), further reducing the symmetry. The space group becomes P 4 (Nr 81), with Cu in (2g), (b) and (c) positions, In in (2g) and (d), Ga in (a) and Se in two distinct (4h) positions. Beyond the Se positions whose all coordinates are not fixed by symmetry, the Cu and In atoms in the (2g) positions have adjustable z coordinates. Obviously, the swap of In↔Ga positions would give rise to the Ga-rich In:Ga=1:3 model supercell.

Figure 2 . 2 :

 22 Figure 2.2: Tetragonal bodycentered supercell used to simulate the In:Ga=1:1 composition in the mixed chalcopyrite-type structure (left panel, the kesterite structure); tetragonal supercell with the In:Ga=3:1 composition (right panel). The (112) planes are marked by color.

2. 3 .

 3 Figure 2.3: BZ and the path selected for energy bands plotting in the (bct) chalcopyritetype structure (left panel) and in the simple tetragonal structure with the same c/a ≈ 2 relation (right panel). The labelling of high-symmetry points is given according to Bradley and Cracknell [1972]. The path chosen is the same in both cases in spite of different labelling.

  [2009], Pohl and Albe [2010] and Oikkonen et al. [2011] used the Perdew-Burke-Ernzerhof (PBE) flavour of the GGA and the Heyd -Scuseria -Ernzerhof (HSE) hybrid XC functional. The values of optimized lattice parameters and the corresponding optical gap, as reported in this works, are reproduced in Table

  the cation and anion sites can be occupied, upwards plane by plane. Reserving capital letters to mark cation positions and lowercase letters for anion positions, moreover retaining the dash to indicate the bonds along [0001], we can write for the wurtzite structure:

Chapter 2 Figure 2 . 4 :

 224 Figure 2.4: Side and top view of the unit cells of β, γ, δ and ε polytypes of GaSe. Large circles: cations, small circles: anions. The stacking notation is explained in the text.

Figure 2 . 5 :

 25 Figure 2.5: A schematic side view of the δ polytype, with displacements of atom layers from highsymmetry positions indicated (strongly exaggerated) by arrows. The labelling of double layers is indicated on the left, and the (x, y) coordinates of corresponding atoms are marked below. The anions of the "Ac" layer are supposed to be in symmetric positions with respect to z=0 / z = 1 2 planes. Cations are shown by small circles and anions by large circles, contrary to the convention of Fig. 2.4.

Figure 2 . 6 :

 26 Figure 2.6: Side view of the hexagonal BZ.

Figure 2 . 7 :

 27 Figure 2.7: Side view of the "fiat" rhombohedral BZ (Figure 1-4 from the PhD thesis of S. Limpijwnnong, Case Western Reserve University, 2000) .Note that the choice of Cartesian unit vectors kx, ky differs by 90°-rotation from that of Kx , Ky used in Sec. 2.8.

2. 5 .

 5 Crystal structures of [In,Ga]Se polytypes column chalcogenides: • GaS crystallizes only in the 2H β [Kuhn et al., 1976]; a=3.587 Å, c=15.492 Å, whence c=7.75 Å. Its internal coordinates as refined by Kuhn et al. [1976]: (

2. 6 .

 6 Crystal structure of In 2 Se 3 tors as layered systems. Thus, Côté et al. [1998] calculated structural and electronic properties of (hypothetical) GaSe nanotubes; Zhirko et al. [2007] discussed perspectives of hydrogen storage in γ-InSe and ε-GaSe; Wang et al. [2009] described highperformance photodetectors whose active element is a InSe nanowire; Late et al. [2012] described a layer transistor based on a single layer of GaS or GaSe; Lei et al. [2014] overview photoconductivity properties of InSe, in view of its use in optoelectronics; Li et al. [2014] described controlled growth of crystalline thin layers of GaSe, potentially useful for electronic and optoelectronic devices.

Figure 2 . 8 :

 28 Figure 2.8: Relation between the wurtzite primitive cell (spanned by black vectors) and the VOSF supercell (spanned by red vectors) in their projection onto the hexagonal plane.

2. 8 .

 8 Appendix 1. BZ of body-centered and simple tetragonal lattices 2.7 Appendix 1. Brillouin zones of body-centered tetragonal and simple tetragonal lattices.Two cases of our interest are the following:1. The bct lattice with c ≈ 2a; translation vectors:a = (-a/2, a/2, c/2 ) ; b = (a/2, -a/2, c/2 ) ; c = (a/2, a/2, -c/2 ) .The reciprocal lattice vectors are:a * bct = (0,2π/a, 2π/c) ; b * bct = (2π/a, 0, 2π/c) ; c * bct = (2π/a, 2π/a, 0) . 2. The simple tetragonal lattice with c ≈ 2a; translation vectors: a = (a, 0, 0) ; b = (0, a, 0) ; c = (0, 0, c) , and the reciprocal lattice vectors a * st = (2π/a, 0, 0) ; b * st = (0, 2π/a, 0) ; c * st = (0, 0, 2π/c) .

Figure 2 . 11 :

 211 Figure 2.11: Reciprocal lattice nodes (blue dots) within the (k x , k z ) plane (above) and within the (k x , k y ) plane (below), along with the corresponding projections of the body-centered tetragonal BZ of the chalcopyrite structure (left panel) and of the simple tetragonal BZ with the same unit cell dimensions (right panel).

Figure 2 . 12 :

 212 Figure 2.12: Projections of the "flat" rhomboherdal BZ with some high-symmetry points indicated; in the rhombohedral setting (left panel); in the hexagonal setting (right panel). See text for details.

aFigure 3 . 2 :

 32 Figure 3.2: Lattice constants, a (red circles) and c/2 (blue triangles) versus composition for Cu(In,Ga)Se 2 . Left panel: WIEN2k calculations with GGA -PBEsol. Experimental values by Hahn et al. [1953] for end systems are indicated by open symbols. Right panel: Fig. 39 from Schnohr [2015]. The data there are taken from Refs. 36 [Suri et al., 1989] and 38 [Durante Rincón et al., 2001]. The dashed and solid lines represent the ternary values and Vegard's law, respectively.

Figure 3

 3 Figure 3.3: Element-specific In-Se (blue triangles), Ga-Se (green circles), and Cu-Se (red squares) bond lengths as a function of composition x for CuIn x Ga 1-x Se 2 . Left panel: WIEN2k calculations with GGA -PBEsol. Right panel: Fig. 41 of Schnohr [2015] depicting the EXAFS measurements. Full and open symbols correspond to powder and polycrystalline samples, respectively; for the meaning of the lines drawn refer to the work cited.

Figure 3 . 4 :

 34 Figure 3.4: Band dispersions as calculated by WIEN2k with PBEsol for CGSe (left panel, showing also the total density of states) and CISe (right panel). Zero energy delimits the occupied bands (shown in blue); the vacant bands are shown in orange.

Figure 3 . 5 :

 35 Figure 3.5: Orbital contributions along band dispersions in CISe as calculated by WIEN2k with PBEsol-GGA. See text for details.

Figure 3 . 6 :

 36 Figure 3.6: Similar to Fig. 3.5 (for the same equilibrium structure), with the mBJ XC potential.

3. 3 .Figure 3 . 7 :Figure 3 . 8 :

 33738 Figure 3.7: Similar to Fig. 3.4 (for the same equilibrium structure), with the mBJ.

Figure 3 . 10 :

 310 Figure 3.10: Left panel: Band gaps in Cu(In,Ga)Se 2 supercells from WIEN2k calculations with the mBJ XC potential (the accurate values are indicated near every data point). Right panel: Fig. 40 from Schnohr [2015], "Bandgap energy versus composition x for the direct semiconductor CuIn x Ga 1-x Se 2 ", citing the ellipsometry results of Durante Rincón et al. [2001].

Fig. 4 .

 4 Fig.4.1 offers a "bird's view" over the performance of different XC schemes, incorporated in WIEN2k and VASP calculations, in what regards the accuracy of predicting the lattice parameters in GaSe and InSe. The polytypes (discussed in Chapter 2) all have similarly small (minimal) lattice parameter in the hexagonal plane, a, but different number of double layers. For the sake of comparison, the c values are "normalized" in the following to the number of structural units per unit cell (two for β and ε, three for γ, four for δ. This will be further on applied also to the discussion of total energies). Moreover it can be mentioned that the WIEN2k calculations for the γ phase have been technically done for the rhombohedral unit cell whereas the VASP calculations used the hexagonal unit cell, tripled with respect to the rhombohedral one. The experimental results are only available for several phases (β,ε-GaSe and β,γ-InSe) studied at different conditions; however they reveal the closeness of a, c parameters throughout the compound given (as well as the closeness of results from different experiments). The scattering of the calculation results, also those obtained within the present work with different XC schemes, occurs at much larger scale.First it should be noted that the GGA-PBE results ("original" PBE, see Chapter 1) are very far off, considerably overestimating both a and c parameters. This can be "repaired" by switching to the PBEsol scheme (applied within the WIEN2k code), or by including dispersive interaction on top of PBE within one or another variation of the Grimme scheme (indicated PBE+D2, PBE+D3 or PBE+D3-BJ) implemented in the VASP code. (Most of the trial calculations with VASP have been done by Michaël Badawi to whom I owe the completeness of these tests). It is obvious that a combination of PBEsol with the Grimme correction (PBEsol+D3, PBEsol-D3-BJ), also tried with VASP, too much overbinds in all cases. As a result of practical observation rather than of a profound analysis, the PBEsol (with WIEN2k) and PBE+D2 (with VASP) are retained for the following detailed inspection of band structures, enery/volume curves, etc.A more attentive look at the results of Fig.4.1 reveals that, whereas the situation with GaSe can be considered as satisfactory (good agreement of PBEsol and PBE-D2 results between themselves and with the experiments available), for InSe it is by far not the case. PBE+D2, while correctly predicting c value, yields an apparently too short (by ∼1%) in-plane lattice constant; PBEsol is performing much better in this sense. The PBEsol, in its turn, has a strange problem with -of all phases -ε-InSe, for which

Figure 4 . 1 :

 41 Figure 4.1: Lattice parameters of GaSe (upper row) and InSe (lower row) polytypes from calculations and experiment. PBE and PBEsol refer to WIEN2k calculations, the cases marked +D2 or +D3 or +D3-BJ -to different implementations in VASP of the Grimme scheme to treat dispersive interaction. The data of earlier calculations are indiscriminately marked by red dots.

Fig. 2

 2 Fig. 2 Schematic band structure of Ga Se.

Figure 4 .

 4 Figure 4.2: Schematic band structure of GaSe: Fig. 2 from Gauthier et al. [1984].

Figure 4 . 3 :

 43 Figure 4.3: Band structures of GaSe polytypes in their respective equilibrium structures, as calculated by WIEN2k with GGA-PBEsol. Occupied (vacant) bands are shown in blue (orange).

Figure 4 . 4 :

 44 Figure 4.4: Band structures of γ-GaSe polytype calculated by WIEN2k with GGA-PBEsol in the rhombohedral setting, along the k-path shown in Fig.2.6 (as for the hexagonal BZ), the same path shifted by [0, 0, ± 2π 3c ], and the superposition of the three band structures into that in the hexagonal setting (right bottom panel, the same as in Fig.4.3). See text for details.

  : WIEN2k(PBEsol) ; VASP (PBE+D2)

Figure 4 . 5 :

 45 Figure 4.5: Comparison of band structures of β-GaSe as calculated by WIEN2k (blue lines) and VASP (red lines). The energy scales are arbitrarily shifted for better visibility. See text for details.

Figure 4 . 6 :

 46 Figure 4.6: Similar to Fig. 4.3 (polytypes in the same relaxed structures), with band structures calculated by WIEN2k with mBJ.

Figure 4 . 7 :

 47 Figure 4.7: Band structures of InSe polytypes in their respective equilibrium structures, as calculated by WIEN2k with mBJ XC potential. Occupied bands are shown in blue, vacant ones -in orange.

4. 4 .

 4 Band structures and Densities of States

Figure 4 . 8 :

 48 Figure 4.8: Band structure with orbital contributions for β-GaSe (upper panel) and InSe (lower panel), as calculated with mBJ.

Figure 4 . 9 :

 49 Figure 4.9: Energy/volume curves as calculated by VASP (left panel) and WIEN2k (right panel)

  of[START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF]:A-C-A• • • B-A-B• • • C-B-C• • • ("Model 1"), or A-B-C• • • B-C-A• • • C-A-B• • • ("Model 2"). Indicating a place of a missing atom in parentheses, this would yield:A-C-A••(C)••B-A-B••(A)••C-B-C••(B)••, or A-B-C••(A)••B-C-A••(B)••C-A-B••(C)••.

Figure 5 . 2 :

 52 Figure 5.2: Band dispersions in three phases of In 2 Se 3 as calculated by WIEN2k with PBEsol.

Figure 5 . 3 :

 53 Figure 5.3: Total (left upper panel) and partial densities of states in LSP1 phase of Ga 2 Se 3 calculated by WIEN2k with PBEsol. See Fig. 5.1 for identification of atoms.

Figure 5 . 5 :

 55 Figure 5.5: Band dispersions in VOSF In 2 Se 3 with s and p orbital contributions from 2coordinated and 3-coordinated Se atoms, as extracted from WIEN2k calculations with PBEsol.

Figure 5 . 6 :

 56 Figure 5.6: Band structures of Ga 2 Se 3 (two left panels) and In 2 Se 3 (two right panels) in the VOSF phase, calculated with PBEsol and mBJ XC potentials.

  in Table2.1. The deviation from the "ideal" value x = 1/4 means that distance of Se to its (two) Cu neighbours is not the same as to its (two) In or Ga neighbours. This distortion is clearly seen in the top view of the (really relaxed) CISe unit cell. The Cu-Se bond is shorter than (In/Ga)-Se, hence x < 1 4 . Mandel et al. [1977] reported the crystal structure of pure CGSe by X-ray diffraction; Paszkowicz et al. [2004] refined the crystal structure of pure CISe. A recent experimental determination of the crystal structure data (by EXAFS) for Cu(In,Ga)Se 2 throughout the whole concentrations range has been recently reported by Schnohr et al. [2012].

	1 4	1 8 ); (x 3 4	1 8 ); ( 1 4	x 7 8 ); ( 3 4 x 7 8 )	Se
	internal coordinate, x

  In have been considered; the corresponding structure considerations are discussed below in Section 2.2. Other structures addressed in this work, the binaries (III)(VI) and (III) 2 (VI) 3 , will be explained in Sections 2.4 and 2.6, respectively. These binary semiconductors belong to the structure phase diagram "bracketing" the (I)(III)(VI) 2 -see, e.g., Figs.3-6 of Pamplin et al. [1979] -or appear in the technical discussions about the sample preparation for photovoltaic devices [Gabor et al., 1996]. Zhang et al. [1998], in their discussion on the defect physics of CISe, invoke the ordered-vacancies phases. In what regards specifically the In-containing phases, these are CuIn 3 Se 5 (1-3-5) and CuIn 5 Se 8 (1-5-8); they seem to be particularly stable and often present in the preparation of Cu-In-Se mixed systems. It was noticed that these stable phases have the generic formula CuIn 2n+1 Se 3n+2 , so that n=0 yields CuInSe

	Chapter 2
	tems with 1:3, 2:2, 3:1 relation of Ga:

2 , n=1 -(1-3-5) and n=2 -(1-5-8). Otherwise we can write generic formula as

1 2 (Cu 2 Se)+(n +

1

2 )(In 2 Se 3

Table 2 . 2 :

 22 Lattice parameters of CISe and CGSe as obtained in earlier calculations within the LDA, in comparison with X-ray diffraction data. a, c: tetragonal lattice parameters; x: the internal coordinate. Crystal structures of [In,Ga]Se polytypes Table 2.3: Results of earlier calculations done on CISe and CGSe with the VASP code, using GGA-PBE and the HSE hybrid functionals. a, c: tetragonal lattice parameters; E g : band gap.

	a ( Å)	c ( Å)	x	Method	Reference
				CuInSe 2
	5.701 11.464 0.2168 VASP (LDA)	Kılıc and Zunger [2003]
	5.733 11.40	0.250 WIEN97 (LDA)	Belhadj et al. [2004]
	5.781 11.609 0.2281 X-ray powder	Zahn and Paufler [1988]
				CuGaSe 2
	5.513 10.941 0.2430 VASP (LDA)	Kılıc and Zunger [2003]
	5.542 10.840 0.260 WIEN97 (LDA)	Belhadj et al. [2004]
	5.5963 11.0036 0.2423 X-ray single crystal Abrahams and Bernstein [1974]

a=5.692 Å, c=11.30 Å end E g =0.38 eV with PBE, and a=5.739 Å, c=11.31 Å end E g =1.73 eV with B3LYP.

Table 3 . 1 :

 31 Lattice parameters of CGSe and CISe according to different experiments.

Table 4 . 1 :

 41 Optimized structure parameters (in Å) of GaSe and InSe polytypes in comparison with experiment: lateral lattice constant a, double layer period c/n and thickness h, vdW gap separating double layers d, cation-anion bond lengths l. For δ phase, average values over a spread of h, d, l present in the structure are given. GaSe-β 3.751 a 8.093 a 4.784 a 3.309 a 2.463 a

	Polytype	a	c/n	h	d	l

Table 4 . 3 :

 43 Band gaps (eV) in GaSe and InSe after GGA-PBEsol and mBJ calculations, in comparison with experimental data.

			Polytypes
	Method	β(2H) γ(3R) ε(2H) δ(4H)
	GaSe PBEsol (Γ-Γ):	0.934	0.924 0.745 0.853
	mBJ (Γ-Γ):	2.092	2.113 1.889 2.010
	mBJ (Γ-M):	1.949	1.963 1.786 1.886
	Exp. a (direct):	2.169	2.120
	Exp. b (direct):		2.020
	Exp. a (indirect): 2.117	2.065
	Exp. b (indirect):		2.010
	InSe PBEsol (Γ-Γ):	0.304	0.240 0.731 0.607
	mBJ (Γ-Γ):	1.232	1.204 1.697 1.493
	Exp. c :		1.29

a Aulich et al. [1969]; b Gauthier et al. [1989]; c Julien and Balkanski [2003].

Table 5 . 2 :

 52 Lattice parameters of Ga 2 Se 3 from experiments.

	a, b, c ( Å)	Identification	Ref.
	5.418	cubic [Hahn α]	Hahn and Klingler [1949]
	6.66, 6.66, 11.65	monoclinic (γ=108.12 • )
		Mikkelsen β	Mikkelsen Jr. [1981]
	5.462	Mikkelsen γ
	cubic 6.661, 11.652, 6.649 monoclinic [Lübbers β] (β=108.94 • ) Lübbers and Leute [1982] Hotje et al. [2005] 5.462 -5.465
	6.645, 11.661, 6.661 monoclinic (β=108.94 • )	Katerynchuk et al. [2014]

There are known cases where the adiabatic approximation cannot be reasonably used, such as dynamics of light atoms, or aspects of photochemistry.

the size of this confinement is controlled by a "cutoff radius" which however cannot be set very small, because this will deteriorate other crucial properties of the pseudopotential. In practice, the markedly "wrong behaviour" of the pseudo-wavefunction may start within ∼10-15% of the interatomic distance.

https://en.wikipedia.org/wiki/Projector_augmented_wave_method

Correspondingly, the "self-designation" of the WIEN2k method is "An Augmented Plane Wave

A scenario of molecular dynamics, in which the atoms follow the forces acting on them but never come to equilibrium, was not applied in the present work.

The γ structure was coded "ABCAB...", the δ one -"ABACA..."

Acknowledgements

, that gave me, during the last year of my thesis, the experience to expand my previous abilities in managing the calculation techniques and to complement the studies done hitherto using WIEN2k.

I gratefully acknowledge the funding sources that made my Ph.D. work possible. I was funded by the Lebanese CNRS fellowship. My work was moreover supported in many possible aspects by the LCP-A2MC Laboratory at the Université de Lorraine and by the Ecole doctorale EDST at the Lebanese University.

List of figures

ε phase includes rotations by 60 • when translating from one double layer to the next one; consequently the "holes" are closed.

Table 2 of [START_REF] Kuhn | Crystal structure and interatomic distances in GaSe[END_REF] summarizes lattice parameters of different polytypes of GaSe, whose structures are shown in Fig. 2. 4. In all of them, a ≃ 3.74 -3.75 Å, and c varies according to the number of planes in the stacking sequence (2 for β and ε; 3 for γ; 4 for δ), so that c, the Z-period per Se-Ga-Ga-Se layer, remains within 7.94 -7.99 Å. Assessing this crystallographic information, one can note that β and ε are the two "simplest" phases, characterized by the shortest possible periodicity along c and

2 ) Ga : 0.180 = Chapter 2 atom in question must be free to move along [START_REF] Gygi | Electronic-structure calculations in adaptive coordinates[END_REF]. Its neighbours in the lower plane (z = 5/6) impose another mirror symmetry and a possibility to displace along [10].

Note that the displacemenmt along the z direction is not fixed either. Taken together, this destroys any symmetry constraints (apart from having 6 partners, related by screw rotation). A practical consequence for choosing calculation input: start from "nominal" 0/ 1 3 / 2 3 positions as in Table 2.5, they are not more "symmetric" than any other, and let the relaxation find the final ones.

In practical calculation setup for either WIEN2k or VASP method, the vacancies do not enter the list of atoms, and any "genuine" atom (2×In + 3×Se) in a given (say z=0) plane makes a separate species, replicated 6 times by screw operations. When passing the coordinates as defined above, the sgroup program in the WIEN2k initialization sequence modifies the input: it shifts the above coordinates by ( 1 3 2 3 0), so that the screw axis passes through (0, 0) -see Sec. 2.9.

The corresponding plane-by-plane distribution of cations (with the space group P 6 1 ) is shown in Fig. 2.10.

The Layer Structure Phase

According to [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF], this phase has wurtzite basal vectors a w 1,2 of Fig. 2.8, with a=4.00 Å and c=28.80 Å. The latter roughly corresponds to 9 wurtzite (double) layers, of which the VOSF phase has 6. Correspondingly, 28.80×2/3=19.20 Å, quite close to 19.38 Å of the VOSF phase but slightly smaller, because now several full layers of In are missing. The 9 double layers go as 3×[-(Se-In-Se-In-Se)--] -see Fig. 10 of [START_REF] Ye | Crystal structures and phase transformation in In 2 Se 3 compound semiconductor[END_REF]; the unit cell has 3 formula units hence 15 atoms.

We can otherwise look at it as an additional anion plane added at the middle of the double-layer block Se-In-In-Se of a InSe structure.

The ambiguity yet remains in the stacking of layers. As there are three layered blocks per unit cell, they can be only stocked as A-B-C. However, within each block the sequence of three Se layers may be either "wurtzite-like", A-B-A, or "zincblende-like", A-B-C. Taken together, this yields two possibilities, shown in Fig. 10

). Indicating a place of a missing atom, this would yield:

Choosing A and B positions as in Fig. 2.8. with C at (0,0), yields the following "ideal" atom positions in the two models. Note that the distance between adjacent cation (or adjacent anion) layers is 1/9, whereas the distance between the cation and anion "on top" is 3/4 of it, i.e., 1/12. 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 ATOM -2: X=0.00000000 Y=0.33333333 Z=0.50000000 MULT= 6 ISPLIT= 8 -2: X=0.66666667 Y=0.00000000 Z=0.66666667 -2: X=0.66666667 Y=0.66666667 Z=0.83333333 -2: X=0.00000000 Y=0.66666667 Z=0.00000000 -2: X=0.33333333 Y=0.00000000 Z=0.16666667 -2: X=0.33333333 Y=0.33333333 Z=0.33333333 In2 NPT= 781 R0=0.00001000 RMT= 2.15 Z: 49.0 LOCAL ROT MATRIX:

1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 ATOM -3: X=0.66666667 Y=0.00000000 Z=0.12500000 MULT= 6 ISPLIT= 8 -3: X=0.66666667 Y=0.66666667 Z=0.29166667 -3: X=0.00000000 Y=0.66666667 Z=0.45833333 -3: X=0.33333333 Y=0.00000000 Z=0.62500000 -3: X=0.33333333 Y=0.33333333 Z=0.79166667 -3: X=0.00000000 Y=0.33333333 Z=0.95833333 Se3 NPT= 781 R0=0.00005000 RMT= 2.18 Z: 34.0 LOCAL ROT MATRIX:

1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 ATOM -4: X=0.66666667 Y=0.00000000 Z=0.45833333 MULT= 6 ISPLIT= 8 -4: X=0.66666667 Y=0.66666667 Z=0.62500000 -4: X=0.00000000 Y=0.66666667 Z=0.79166666 -4: X=0.33333333 Y=0.00000000 Z=0.95833333

Chapter 3 gration of the wavefunctions within the (non-overlapping) atomic radii, therefore the contribution from the interstitial region is missing. Moreover, the contributions in some other l channels than those shown in Fig. 3.8 (e.g., Cu-p, Se-d etc.) are not negligible. In order to check how the total DOS sums up from partial ones one should remember to count the Se contribution twice, moreover multiply the sum by two (because of two formula units in the primitive cell) and once more by two (because of spin degeneracy).

Band structures in mixed systems; band gaps varying with concentration

Fig. 3.9 depicts the energy bands for superlattices which imitate intermediate concentrations of Ga and In. As discussed in Chapter 2, the 50:50 concentration is modeled by letting the two III-column cation positions in the chalcopyrite unit cell to be occupied one by Ga and the other by In. The Brillouin zone and the number of bands are not changed as compared to the case of pristine compounds, however, some degeneracies in the band dispersions are lifted. This is best seen by comparing the bands along the X-P path or in Z, where many degeneracies occured in Fig. 3.7 which are lifted as the In/Ga distinction happens between the two sites. For one thing, the Ga4s and In5s Chapter 4 phases reveals a lot of degeneracies in A, which are lifted along the (A-Γ) line. This is the most clearly seen for the isolated bands near -5 --7 eV, and also in the highest occupied bands. If such "simple" splitting occurs in β, ε systems, in γ we observe a more sophisticated "back and forth" splitting, due to tripling the structure unit, and in δ -the "reflected" splitting, as a band originally degenerate in Γ splits towards A, and each of its parts splits yet again on returning from A to Γ. In the lack of interaction between double layers such folding would be just a geometry exercise; in real systems, however, the upfolded bands are slightly modified by interlayer interaction, and namely there modifications give rise to tiny preferences in the band energy of one or the other phase, to be discussed in the last section of this chapter.

Another consequence of the foldings of bands is a possible modification of the band gap character, as the positions of the bands' minima and maxima may be affected by replicating the unit cell. For example, the band gap which occurs at Γ in the β and ε phases shifts to A in the γ phase, as an even number of bands is affected by the "back and forth" splitting. One can imagine that in some other possible polytypes, the nature of the band gap might be changed from direct to indirect.

Specifically for the γ phase, the stockpiling of bands can be the most naturally analyzed. The calculation for this phase has been done in WIEN2k as for the rhombohedral structure, with one double layer (four atoms) per unit cell. The counting of the bands in this case, following from the counting of electrons included in the valence states (above the atomic cores), is the following: Se 3d 10 4s 2 4p 4 → 16 electrons + Ga 3d 10 4s 2 4p 1 → 13 electrons = 29 × 2 (formula units per unit cell) = 58 electrons → 29 bands (of which the deepest ten, due to Se3d states, are not covered by the band structure plots shown; the next ten, due to Ga3d, form a tight bunch near -15 eV, and the remaining nine can be easily counted). This relatively "sparse" band structure is then traced along three different k paths in Fig. 4.4, according to how our "standard" hexagonal path is mapped onto the (three times larger) Brillouin zone of the rhombohedral lattice. In fact this amounts to simply shifting the path from its original position by ±[0 0, 2π/(3c)]. The three resulting "rhombohedral" band structures can be then simply superposed, as is shown in Fig. 4.4, where the color code of each contributing figure is kept different, for clarity. The resulting band structure plot is exactly the same as would be obtained from performing a calculation on an equivalent, three times larger hexagonal unit cell. The corresponding relations between unit cells and Brillouin zones are explained in the Appendix of Chapter 2.

For making sure that the band structures delivered by WIEN2k and VASP are consistent, a direct comparison of the results of two calculations, for the case of β-GaSe, is shown in Fig. 4.5. The underlying crystal structures (for which the bands were calculated) stem from the relaxation with PBEsol (in case of WIEN2k) and PBE+D2 (in case of VASP) schemes. As can be seen from Fig. 4.1 and Table 4.1, the resulting structures are reasonably close. The band were calculated using, correspondingly, the The second set of curves in the left panel of Figure 4.9, marked "c/a relaxed", which share the minimum points with their "c/a fixed" counterparts but pass much more flat, represent the case where the structure was fully relaxed (including also the internal coordinates) at each trial value of volume.

These curves could be used for extracting the bulk modulus (not discussed here because of limited usefulness of this single parameter for the strongly anisotropic systems under study). Corresponding calculations have not been done with WIEN2k, because of too large amount of relaxations to be done effectively "by hand".

Similar calculations (identification of the fully relaxed structure for each phase, followed by "rigid" and "c/a relaxed" volume scans) have also been done for InSe; the results (from VASP calculations only) are shown in Fig. 4.4. Differently from the case of GaSe, one notes a small spread of minimum positions -from ∼117.5 Å3 ) (for γ and ε, which are almost degenerate and come about as ground states) to ∼118.5 Å3 ) for the β phase which is now characterized by the highest-energy curve, opposite to its attribution in GaSe. One notes however that at small negative pressure (i.e., looking at the range of volumes larger than 120 Å3 ), the β phase would tend to dominate, and, in general, the whole system of E(V ) curves very faithfully recovers the order of phases established for GaSe. We note that the predominance of the γ phase for InSe was long ago argued for by [START_REF] Likforman | Diagramme de phases du système indiumsélénium[END_REF].

Chapter 5

Results for Ga 2 Se 3 and In 2 Se 3

Introduction about crystal structures

An attempt to address the In 2 Se 3 and Ga 2 Se 3 systems faces a problem of certain ambiguity in crystal structures definition. The statement by [START_REF] Popović | X-ray diffraction study and semiconducting properties of the system Ga 2 Se 3 -In 2 Se 3[END_REF], The published data concerning In,Se, and its phase transitions are rather confusing and even conflicting in many details, is still valid nowadays. In principle, the (III) 2 (VI) 3 stoichiometry restores the formal matching of cation / anion valences, so that the structure can be apparently arranged without "neutralizing" an excess in cations by forming cation-cation bonds, using just a network of "regular" cation-anion bonds, with, presumably, tetrahedral local environments. On the other hand, the 2:3 stoichiometry can come together with tetrahedral connectivity only if 1/3 of cation vacancies are somehow included into the lattice, according to the formula (III) 2 (VI) 3 . What is known with relative safety is that the hexagonal (or rhombohedral) setting is predominant, and the sequence of Se-In-Se-In-Se-layers is an important structure element; however, the stoichiometry can be established by placing the cation vacancies in the layers, instead of organizing a full missing layer.

Table 5.1 offers some limited overview of the experimental data concerning In 2 Se 3 . In total, we assume that the VOSF and LSP models introduced by Ye et al. [1998] and discussed in detail in Sec. 2.6 are plausible from the side of experiment, and make a valid basis for theory simulations.

For Ga 2 Se 3 , there are largely other structures than hexagonal that are brought into consideration -see Table 5.2. Most of the later work refer to phases which can be identified with those indicated as "Hahn α" and "Lübbers β" structures, both being stable. The most thorough analysis of available crystal structure data, among relatively recent works, has been offered by [START_REF] Hotje | Chemischer Transport fester Lösungen[END_REF] and [START_REF] Huang | First-principles study of the structure and band structure of Ga 2 Se 3[END_REF]. These latter authors essentially reported first-principles studies of the electronic structure, optimizing also the lattice parameters, for monoclinic and orthorhombic phases. Much Chapter 5 nates). The resulting parameters are given in Table 5.3. The lattice parameters of LSP and VOSF phases are very different, but both are derived from the wurtzite lattice (by introducing vacancies either in "scattered" form, as in VOSF, or as missing cation planes in LSP). It is convenient to refer to underlying wurtzite parameters, which can be recovered as a/ √ 3, c/6 for VOSF and a, c/9 for LSP.

Relaxed crystal structures from calculation

The Wyckoff positions of the relaxed structures, along with corresponding band lengths, are listed in Tables 5.4, 5.5 for the VOSF phase and in Tables 5.6, 5.7 for the layered phases. It is obvious that the bond connecting the "lonely" Se atom (2coordinated in the VOSF phase, or "top Se" in layered phases) is considerably shortened.

As follows from the discussion in Sec. 5.1, the comparison of calculated lattice parameters with experiment is not possible in case of Ga 2 Se 3 . However, an interesting comparison can be made with the work by [START_REF] Popović | X-ray diffraction study and semiconducting properties of the system Ga 2 Se 3 -In 2 Se 3[END_REF], who studied the variation of structures in (In,Ga) 2 Se 3 solid solutions. According to their data, the mixed crystals remain hexagonal starting from the In-rich side all the way through up to 1:2 molar ratio of In:Ga, whereby the hexagonal a parameter gradually shrinks from 7.11 Å in In 2 Se 3 to 6.82 Å in (In 1/3 Ga 2/3 ) 2 Se 3 . At the same time, the c parameter starts from 19.33 Å at the In-reach side, gradually drives through the minimum c = 18.85 Å at In:Ga=1:1 concentration, and recovers to c = 19.30 Å at the Ga concentration reaches 2/3 (and the hexagonal structure ceases to exist). It is interesting to note that the corresponding lattice parameter values from Table 5.3 are quite close to these numbers, including the fact that the c parameters do not differ much between Ga 2 Se 3 and In 2 Se 3 (within ∼1.5%), whereas the a parameters do (by ∼6%). It would be interesting and not particularly difficult to estimate, from additional calculations for the 1:1 mixed system, whether the c parameter will be indeed considerably reduced, as compared to end systems, and to try to understand the driving forces of such effect.

As could be expected, the in-plane size of VOSF phase, that contains vacancies in every layer, is more tight. Interestingly, the "wurtzite c" parameter of LSP phases is (as intuitively expected) increased in Ga 2 Se 3 but (a bit counter intuitively) expanded in In 2 Se 3 , with respect to that of the VOSF phase. This is consistent with the fact that the LSP structures are hold together by interaction between triple layers across the vdW gap, whereas the VOSF structures are knitted by "conventional" covalent interactions in all directions. The dispersive interactions, a priori not included in a conventional The corresponding band structures, in comparison according to PBEsol and mBJ calculations, are shown in Fig. 5.6. The band structure of In 2 Se 3 with PBEsol was already shown above in Fig. 5.2 but is included in Fig. 5.6 for the convenience of comparison. One can see that the modifications of the band structure are more complex than just increasing the band gap. Namely, the (upper) valence band (dominated by Se4p and Ga4p/In5p states, see Fig. 5.4) gets "compressed" by ∼0.5 eV, "pushed" from below by the nearly rigid bunch of Se4p-Ga4s/In5s states. The Se4s-dominated bands situated at about 13 eV below the valence band top remain essentially rigid on passing from PBEsol to mBJ; they only undergo a slight upward shift. The Ga3d / In4d-related band at ∼15 eV, on the contrary, drifts downwards and becomes more narorow in mBJ. The same effect of the Se4s band shifting upwards against the In4d band shifting downwards under the effect of the mBJ XC potential can be seen in Fig. 3.5 vs Fig. 3.6 in Chapter 3.

The conduction band, on being generally shifted upwards under the effect of mBJ, becomes more flat. Whereas the top of the valence band stays invariably in Γ, a local minimum in the conduction band develops at M and becomes competitive with that in Index of abbreviations