Appendix A Ab initio Computational Methods The objective of ab initio approaches based on applied quantum theory is to calculate stationary states for electrons in the electrostatic field of nuclei,i.e., the electronic structure (ES).The energy of this ground state can then serve as a potential energy for displacements of nuclei.From the point of view of ideal-strength (IS)calculations,the total energy of the system is the most important output of ab initio methods. The main advantage of ab initio methods is their independence of exper- imental data.Unlike in the case of empirical and semi-empirical methods, there is no need for calibration parameters.Thus,they can also be used for calculation of some structural and mechanical characteristics of hypothetical systems(prediction of properties of materials that have not yet been devel- oped)or study of materials behaviour under large deformations(far from the equilibrium state)that can give a better understanding of micromechanisms of materials failure. First attempts to develop applicable theories were made in the late 1920s [421,422,a few years after the foundations of modern quantum the- ory were laid(derivation of the Schrodinger equation).A very successful step forward was the Hartree-Fock method [422,423].This method yields very accurate bond lengths in molecules.On the other hand,the binding energies are generally not in good agreement with experimentally obtained energies. Moreover,for solids,the Hartree-Fock method has problems with a descrip- tion of band structures.The density functional theory (DFT)37,38]was invented to include correlation effects without using the very costly wave function methods.All the methods used within this book are based on the DFT. In DFT the energy is not obtained as eigenvalues of a wave function,but rather as a functional of the electron density.The complex problem of many interacting electrons is transformed into a much simpler study of single elec- tron interactions with an effective potential Uef created by other electrons and all nuclei.This is expressed by the Kohn-Sham equation (one-electron Schrodinger equation) 249
Appendix A Ab initio Computational Methods The objective of ab initio approaches based on applied quantum theory is to calculate stationary states for electrons in the electrostatic field of nuclei, i.e., the electronic structure (ES). The energy of this ground state can then serve as a potential energy for displacements of nuclei. From the point of view of ideal-strength (IS) calculations, the total energy of the system is the most important output of ab initio methods. The main advantage of ab initio methods is their independence of experimental data. Unlike in the case of empirical and semi-empirical methods, there is no need for calibration parameters. Thus, they can also be used for calculation of some structural and mechanical characteristics of hypothetical systems (prediction of properties of materials that have not yet been developed) or study of materials behaviour under large deformations (far from the equilibrium state) that can give a better understanding of micromechanisms of materials failure. First attempts to develop applicable theories were made in the late 1920s [421, 422], a few years after the foundations of modern quantum theory were laid (derivation of the Schr¨odinger equation). A very successful step forward was the Hartree–Fock method [422, 423]. This method yields very accurate bond lengths in molecules. On the other hand, the binding energies are generally not in good agreement with experimentally obtained energies. Moreover, for solids, the Hartree–Fock method has problems with a description of band structures. The density functional theory (DFT) [37, 38] was invented to include correlation effects without using the very costly wave function methods. All the methods used within this book are based on the DFT. In DFT the energy is not obtained as eigenvalues of a wave function, but rather as a functional of the electron density. The complex problem of many interacting electrons is transformed into a much simpler study of single electron interactions with an effective potential Ueff created by other electrons and all nuclei. This is expressed by the Kohn–Sham equation (one-electron Schr¨odinger equation) 249
250 A Ab initio Computational Methods (-△+Uef(r)-e)(r)=0, where si represents the one-electron energies and i are the one-electron wave functions. The wave functions are then occupied in accordance with the Pauli prin- ciple and a new field is obtained by solving the Poisson equation for point nuclei shielded by the electronic density p(r)=>lv:(r)2. i.occ In the case of periodic crystalline materials,the one-electron wave func- tions are expanded into appropriate basis sets and satisfy the Bloch theorem. The Kohn-Sham equation is solved iteratively until the solution becomes self-consistent.This means that the electron density,determined from the effective one-electron potential,must generate the same effective potential (which is again a functional of the electron density).The self-consistent cycle usually starts with a guess of the effective potential(superposition of atomic- like potentials)and then the input and output potentials are appropriately mixed before starting a new iteration.The quality and speed of the conver- gence of such calculations is related not only to the choice of a suitable basis but also to the sophistication of the iterative process.The necessary correc- tions for exchange and correlation are also included in the effective potential Uef.This seems to be the crucial point of ab initio calculations because the exchange-correlation(XC)functional is not known exactly and must be ap- proximated.The first (and the simplest)attempt to build an approximation of the XC energy functional in the DFT is the local density approximation (LDA)[424].The LDA is local in the sense that the electron exchange and correlation energy at any point in space is a functional of the electron density at that point only.As a consequence of this,LDA fails in situations where the density undergoes rapid changes (molecules).An improvement to this can be made by considering the gradient of the electron density.The density gradient corrections are implemented in the so-called Generalized Gradient Approximation(GGA).While there is only one LDA there are several differ- ent parameterizations of the GGA [425-427]. Nevertheless,it is the choice of the basis wave functions that makes the main difference between various methods used for the ES calculations.The better we choose them (according to the character of the problem),the smaller number of them is needed for the description of one-electron wave functions.Commonly used bases are augmented (APW)and orthogonalized (OPW)plane waves,linear muffin-tin orbitals (LMTO),linear combination of atomic orbitals (LCAO),of Gaussian orbitals (LCGO)and linear aug- mented Slater-type orbitals (LASTO),augmented spherical waves (ASW), etc.The Korringa-Kohn-Rostoker(KKR)method proceeds by the use of the Green function of the Kohn-Sham equation and is also called Green's func-
250 A Ab initio Computational Methods (−Δ + Ueff (r) − εi)ψi(r)=0, where εi represents the one-electron energies and ψi are the one-electron wave functions. The wave functions are then occupied in accordance with the Pauli principle and a new field is obtained by solving the Poisson equation for point nuclei shielded by the electronic density ρ(r) = i,occ |ψi(r)| 2 . In the case of periodic crystalline materials, the one-electron wave functions are expanded into appropriate basis sets and satisfy the Bloch theorem. The Kohn–Sham equation is solved iteratively until the solution becomes self-consistent. This means that the electron density, determined from the effective one-electron potential, must generate the same effective potential (which is again a functional of the electron density). The self-consistent cycle usually starts with a guess of the effective potential (superposition of atomiclike potentials) and then the input and output potentials are appropriately mixed before starting a new iteration. The quality and speed of the convergence of such calculations is related not only to the choice of a suitable basis, but also to the sophistication of the iterative process. The necessary corrections for exchange and correlation are also included in the effective potential Ueff . This seems to be the crucial point of ab initio calculations because the exchange-correlation (XC) functional is not known exactly and must be approximated. The first (and the simplest) attempt to build an approximation of the XC energy functional in the DFT is the local density approximation (LDA) [424]. The LDA is local in the sense that the electron exchange and correlation energy at any point in space is a functional of the electron density at that point only. As a consequence of this, LDA fails in situations where the density undergoes rapid changes (molecules). An improvement to this can be made by considering the gradient of the electron density. The density gradient corrections are implemented in the so-called Generalized Gradient Approximation (GGA). While there is only one LDA there are several different parameterizations of the GGA [425–427]. Nevertheless, it is the choice of the basis wave functions that makes the main difference between various methods used for the ES calculations. The better we choose them (according to the character of the problem), the smaller number of them is needed for the description of one-electron wave functions. Commonly used bases are augmented (APW) and orthogonalized (OPW) plane waves, linear muffin-tin orbitals (LMTO), linear combination of atomic orbitals (LCAO), of Gaussian orbitals (LCGO) and linear augmented Slater-type orbitals (LASTO), augmented spherical waves (ASW), etc. The Korringa–Kohn–Rostoker (KKR) method proceeds by the use of the Green function of the Kohn–Sham equation and is also called Green’s func-
A.1 TB-LMTO-ASA Code 251 tion (GF)method.The pseudopotential approach applied mostly to solids containing no d-or f-electrons is also widely used.A detailed description of these methods may be found in many books and articles,e.g.,in 86,88,89. The atomic configurations corresponding to deformed structures usually have lower symmetries and,at the strength limit,they are very far from the lowest- energy equilibrium state.Therefore,to obtain reliable structural energy dif ferences,the full-potential methods (i.e.,without any shape approximation of the crystal potential and the electronic charge density)have to be utilized in such studies.At present,several codes are available,e.g.,WIEN,VASP, FHI,FLEUR,FPLO,FPLMTO,ABINIT,SIESTA,etc. All the ab initio calculations of IS presented in this book were performed by using the following three computational codes:LMTO-ASA,WIEN and VASP. A.1 TB-LMTO-ASA Code Linearized ab initio methods have been successfully utilized for solving many problems in materials science [428,429].One of the most effective approaches for early first principles calculations was the LMTO formalism which has been continuously developed since 1980 [428-430.This method is very appropriate for self-consistent calculations. In the LMTO-ASA code [431],the crystal potential U is approximated by a muffin-tin(MT)shape potential which is composed of a set of spherically symmetric potentials inside slightly overlapping spheres around individual nuclei and a constant potential in the interstitial region outside the spheres (Figure A.1).Atomic-like orbitals derived for the MT potential constitute a suitable basis set.The tight binding approximation,which is also imple- mented into the code [432],assumes that the full Hamiltonian may be ap- proximated by that of an isolated atom centred at each lattice point.The atomic orbitals (eigenfunctions of the single-atom Hamiltonian)are assumed to be negligible at distances exceeding the lattice constant. In all presented calculations,the LMTO method is used in the framework of an atomic sphere approximation(ASA)which is particularly suitable for closely packed structures like fcc,hcp or bcc [428].The size of a spherically symmetric potential is assumed to be equal to that of the Wigner-Seitz cell (Figure A.1(b)).This suppresses the interstitial region and neglects the ki- netic energy of related free electrons.Owing to the necessary space-filling condition,ASA represents a physically plausible model only for a description of an infinite system of atomic spheres
A.1 TB-LMTO-ASA Code 251 tion (GF) method. The pseudopotential approach applied mostly to solids containing no d- or f-electrons is also widely used. A detailed description of these methods may be found in many books and articles, e.g., in [86, 88, 89]. The atomic configurations corresponding to deformed structures usually have lower symmetries and, at the strength limit, they are very far from the lowestenergy equilibrium state. Therefore, to obtain reliable structural energy differences, the full-potential methods (i.e., without any shape approximation of the crystal potential and the electronic charge density) have to be utilized in such studies. At present, several codes are available, e.g., WIEN, VASP, FHI, FLEUR, FPLO, FPLMTO, ABINIT, SIESTA, etc. All the ab initio calculations of IS presented in this book were performed by using the following three computational codes: LMTO-ASA, WIEN and VASP. A.1 TB-LMTO-ASA Code Linearized ab initio methods have been successfully utilized for solving many problems in materials science [428,429]. One of the most effective approaches for early first principles calculations was the LMTO formalism which has been continuously developed since 1980 [428–430]. This method is very appropriate for self-consistent calculations. In the LMTO-ASA code [431], the crystal potential U is approximated by a muffin-tin (MT) shape potential which is composed of a set of spherically symmetric potentials inside slightly overlapping spheres around individual nuclei and a constant potential in the interstitial region outside the spheres (Figure A.1). Atomic-like orbitals derived for the MT potential constitute a suitable basis set. The tight binding approximation, which is also implemented into the code [432], assumes that the full Hamiltonian may be approximated by that of an isolated atom centred at each lattice point. The atomic orbitals (eigenfunctions of the single-atom Hamiltonian) are assumed to be negligible at distances exceeding the lattice constant. In all presented calculations, the LMTO method is used in the framework of an atomic sphere approximation (ASA) which is particularly suitable for closely packed structures like fcc, hcp or bcc [428]. The size of a spherically symmetric potential is assumed to be equal to that of the Wigner–Seitz cell (Figure A.1(b)). This suppresses the interstitial region and neglects the kinetic energy of related free electrons. Owing to the necessary space-filling condition, ASA represents a physically plausible model only for a description of an infinite system of atomic spheres
252 A Ab initio Computational Methods 5 0 3 3 1012345 G- [001] [110] -7012 [001) 34 [110] (a) (b) Figure A.1 Muffin-tin shape potential in(110)plane of a general bcc crystal of the lattice constant a =3au with spherical radii:(a)rMT 1au,and (b)rMT =rws 1.48au A.2 Wien 95-w2k Codes The program package WIEN [433]does not use any shape approximation to the potential.The crystal environment is divided into a region of non- overlapping atomic spheres(centred at individual atomic sites)and an inter- stitial region as can be seen in Figure A.2.In order to describe ES reliably and effectively,two different basis sets are employed:a product of linear com- bination of radial functions and spherical harmonics is used inside the spheres whereas the wave functions in the interstitial region are expanded into a linear combination of plane waves.The solution in both regions must be continuous at the spherical boundaries.Each basis function is then defined as a plane- wave in the interstitial region connected smoothly to a linear combination of atomic-like functions in the spheres,thus providing an efficient represen- tation throughout the space.A similar representation is used for potentials and charge densities.The method is called the linear augmented plane wave (LAPW)method [428,434]. S Figure A.2 Illustration of a crystal model-three atomic spheres (S1-S3)with (embedded in the interstitial region I with
252 A Ab initio Computational Methods (a) (b) Figure A.1 Muffin-tin shape potential in (110) plane of a general bcc crystal of the lattice constant a = 3 au with spherical radii: (a) rMT = 1 au, and (b) rMT = rW S = 1.48 au A.2 Wien 95 – w2k Codes The program package WIEN [433] does not use any shape approximation to the potential. The crystal environment is divided into a region of nonoverlapping atomic spheres (centred at individual atomic sites) and an interstitial region as can be seen in Figure A.2. In order to describe ES reliably and effectively, two different basis sets are employed: a product of linear combination of radial functions and spherical harmonics is used inside the spheres whereas the wave functions in the interstitial region are expanded into a linear combination of plane waves. The solution in both regions must be continuous at the spherical boundaries. Each basis function is then defined as a planewave in the interstitial region connected smoothly to a linear combination of atomic-like functions in the spheres, thus providing an efficient representation throughout the space. A similar representation is used for potentials and charge densities. The method is called the linear augmented plane wave (LAPW) method [428, 434]. S1 S2 S3 I Figure A.2 Illustration of a crystal model – three atomic spheres (S1 – S3) with potential US(r) = 2 lm Ulm(r)Ylm(ˆr) embedded in the interstitial region I with UI (r) = 2 K UK(r)eiKr
A.3 VASP Code 253 In order to increase the flexibility of the basis (to improve upon the lin- earization of wave functions)and to make possible a consistent treatment of semicore and valence states in one energy window(to ensure orthogonality) additional basis functions can be added.They are called local orbitals [435 and consist of a linear combination of two radial functions at two different energies (e.g.,at the 3s and 4s energy)and one energy derivative (at one of these energies).The local orbitals are normalized and have zero value and slope at the spherical boundaries. A.3 VASP Code Another way of avoiding a problem with plane wave basis set in the vicinity of atomic nuclei,where the number of plane waves would exceed any practical limits(perhaps except for H or Li),is to substitute the exact potential by a pseudopotential. LU cut Ψ al-electron cut pseudo U- Figure A.3 Comparison of a wavefunction in the Coulomb potential of the nucleus (dashed lines)to that in the pseudopotential (solid lines) Construction of the pseudo-wavefunctions is schematically described in Figure A.3.The Coulomb potential of the nucleus and the corresponding wavefunction are represented by dashed lines.The solid line displays the pseudopotential and the pseudo-wavefunction.The real and the pseudo-
A.3 VASP Code 253 In order to increase the flexibility of the basis (to improve upon the linearization of wave functions) and to make possible a consistent treatment of semicore and valence states in one energy window (to ensure orthogonality) additional basis functions can be added. They are called local orbitals [435] and consist of a linear combination of two radial functions at two different energies (e.g., at the 3s and 4s energy) and one energy derivative (at one of these energies). The local orbitals are normalized and have zero value and slope at the spherical boundaries. A.3 VASP Code Another way of avoiding a problem with plane wave basis set in the vicinity of atomic nuclei, where the number of plane waves would exceed any practical limits (perhaps except for H or Li), is to substitute the exact potential by a pseudopotential. Upseudo r rcut pseudo r rcut U 0 0 Figure A.3 Comparison of a wavefunction in the Coulomb potential of the nucleus (dashed lines) to that in the pseudopotential (solid lines) Construction of the pseudo-wavefunctions is schematically described in Figure A.3. The Coulomb potential of the nucleus and the corresponding wavefunction are represented by dashed lines. The solid line displays the pseudopotential and the pseudo-wavefunction. The real and the pseudo-