計算聲子譜以及熱力學(xué)性質(zhì)的綜述

1、#Baroni, Giannozzi, IsaevReviews in Mineralogy & GeochemistryVol. 71 pp. XXX-XXX. 2009Copyright © Mineralogical Society of AmericaThermal Properties of Materials fromAb Initio Quasi-Harmonic PhononsStefano BaroniScuola Internazionale Superiore di Studi Avanzativia Beirut 2-434151 Trieste, I

2、talybaroniissa.itPaolo GiannozziDipartimento di Fisica, Universita di Udine,Udine CNR DEMOCRITOS National Simulation CenterTriestepaolo. giannozzi iniud. itEyvaz IsaevDepartment of Physics, Chemistry and Biology (IFM)Linkoping UniversitySwedenTheoretical Physics DepartmentMoscow State Institute of S

3、teel and AlloysRussiaeyazjsaev(INTRODUCTIONComputer simulations allow for the investigation of many materials propenies and processes that are not easily accessible in the laboratory This is particularly true in the Earth sciences, where the relevant pressures and temperatures may be so extreme that

4、 no experimental techniques can operate at those conditions Computer modeling is often the only source of information on the properties of materials that, combined with indirect evidence (such as seismic data), allows one to discriminate among competing planetary models. Many computer simulations ar

5、e performed using effective inter-atomic potentials tailored to reproduce some experimentally observed properties of the materials being investigated. The remoteness of the physically interesting conditions from those achievable in the laboratory, as well as the huge variety of different atomic coor

6、dination and local chemical state occurring in the Earth interior, make the dependability of semi-empirical potentials questionable. First-principles techniques based on density-functional theory (DFT) (Hohenberg and Kohn 1964; Kohn and Sham 1965) are much more predictive, not being biased by any pr

7、ior experimental input, and have demonstrated a considerable accuracy in a wide class of materials and variety of external conditions. The importance of thermal effects in the range of phenomena interesting to the Earth sciences makes a proper account of atomic motion essential. Traditionally, this

8、is achieved using molecular dynamics techniques which have been successfully combined with DFT in the first-principles molecular dynamics technique of Car and Parrinello (1985). Well below the melting temperature, the numerical efficiency of molecular dynamics is limited by the lack of ergodicity, w

9、hich would require1529-6466/09/0071 -OOO3$O5.00DOI: 10.2138/rmg.2OO9.71.3 long simulation times, and by the importance of long-wavelength collective motions (phonons), which would require large simulation cells Both difficulties are successfully dealt with in the quasi-harmonic approximation (QHA) w

10、here the thermal properties of solid materials are traced back to those of a system of non-interacting phonons (whose frequencies are however allowed to depend on volume or on other thermodynamic constraints). An additional advantage of the QHA is that it accounts for quantummechanical zero-point ef

11、fects, which would not be accessible to molecular dynamics with classical nuclear motion. The availability of suitable techniques to calculate the vibrational properties of extended materials using a combination of DFT and linear-response techniques (resulting in the so-called density-functional per

12、turbation theory, DFPT; Baroni et al. 1987, 2001) makes it possible to combine the QHA with DFT. The resulting simulation methodology has shown to be remarkably accurate in a wide temperature range, extending up to not very far from the melting line and has been applied to a wide variety of systems,

13、 including many which are relevant to the Earth sciences This paper gives a short overview of the calculation of thermal properties of materials in the framework of the QHA, using DFT. The paper is organized as follows: in the "Thermal Properties and the QuasiHarmonic Approximation section, we

14、introduce some of the thermal properties of interest and describe how they can be calculated in the framework of the QHA; in the "Ab Initio Phonons" section, we describe the DFPT approach to lattice dynamics; in the "Computer Codes" section, we briefly introduce some of the compu

15、ter codes that can be used to perform this task; in the "Applications5' section, we review some of the application of the first-principles QHA to the study of the thermal properties of materials; finally, the last section contains our conclusions.THERMAL PROPERTIES AND THE QUASI-HARMONIC AP

16、PROXIMATIONThe low-temperature specific heat of solids is experimentally found to vanish as the cube of the temperature, with a cubic coefficient that is system-specific (Kittel 1996; Wallace 1998). This is contrary to the predictions of classical statistical mechanics, according to which the heat c

17、apacity of a system of harmonic oscillators does not depend on temperature, nor on its spectrum. One of the landmarks of modem solid-state physics, that greatly contributed to the establishment of our present quantum-mechanical picture of matter, is the Debye model for the heat capacity of solids. T

18、his model naturally explains the low-temperature specific heat of solids in terms of the (quantum) statistical mechanics of an ensemble of harmonic oscillators, which can in turn be pictorially described as a gas of non-interacting quasi-particles obeying the Bose-Einstein statistics (phonons).The i

19、nternal energy of a single harmonic oscillator of angular frequency © in thermal equilibrium at temperature T. is:<E>=where kB is the Boltzmann constant. By differentiating with respect to temperature the sum over all the possible values of the phonon momentum in the Brillouin zone (BZ) o

20、f Equation (1), the constant-volume specific heat of a crystal reads:CV(T) = Mq,v)nz(q.v)V qvwhere co(q,v) is the frequency of the v-th mode (phonon) at point q in the BZ, zf(q,v)= (3 / 3T)e(*<0<,v)/*>r) l_l, and the sum is extended to the first BZ. By assuming that there are three degenera

21、te modes at each point of the BZ, each one with frequency GXq,v) = c|q|, c being the sound velocity, and converting the sum in Equation (2) into an integral, the resulting expression for the heat capacity, valid in the low-temperature limit, reads:Baroni etal indd 212/9/2009 2 15 59 PMThermal Proper

22、ties from Ab Initio Quasi-Harnionic Phonons5where Q is the volume of the crystal unit cell and 0D = (2兀力 / 心)c(3 / 47cQ),/3 is the so-called Debye temperature In the BomOppenheimer approximation (Bom and Oppenheimer 1927), the vibrational properties of molecules and solids are determined by their el

23、ectronic structure through the dependence of the ground-state energy on the coordinates of the atomic nuclei (Martin 2004). At low temperature the amplitudes of atomic vibrations are much smaller than interatomic distances, and one can assume that the dependence of the ground-state energy on the dev

24、iation from equilibrium of the atomic positions is quadratic. In this, so called harmonic, approximation (HA) energy differences can be calculated from electronic-structure theory using static response functions (DeCicco and Johnson 1969; Pick et al. 1970) or perturbation theory (Baroni et al. 1987,

25、 2001) (see the next section).In the HA, vibrational frequencies do not depend on interatomic distances, so that the vibrational contribution to the crystal internal energy does not depend on volume. As a consequence, constant-pressure and constant-volume specific heats coincide in this approximatio

26、n, and the equilibrium volume of a crystal does not depend on temperature Other shortcomings of the HA include its prediction of an infinite thermal conductivity, infinite phonon lifetimes, and the independence of vibrational spectra (as well as related properties: elastic constants, sound velocitie

27、s etc.) on temperature, to name but a few. A proper account of anharmonic effects on the static and dynamical properties of materials would require the calculation of phonon-phonon interaction coefficients for all modes in the BZ Although the leading terms of such interactions can be computed even f

28、rom first principles (Baroni and Debemardi 1994; Debemardi et al. 1995)and the resulting vibrational linewidths have in fact been evaluated in some cases (Debemardi et al. 1995; Lazzeri et al. 2003; Bonini et al. 2007)the extensive sampling of the phonon-phonon interactions over the BZ required for

29、free-energy evaluations remains a daunting task. The simplest generalization of the HA, which corrects for most of the above mentioned deficiencies, while not requiring any explicit calculation of anharmonic interaction coefficients, is the QHAIn the QI I A, the crystal free energy is assumed to be

30、determined by the vibrational spectrum via the standard harmonic expression:F(X,T) = UQ(X) +| X) + kBT 1 一/護(hù) 1/ qvqvV丿where X indicates any global static constraint upon which vibrational frequencies may depend (most commonly just the volume V, but X may also include anisotropic components of the st

31、rain tensor, some externally applied fields, the internal distortions of the crystal unit cell, or other thermodynamic constraints that may be applied to the system), and U0(X) is the zerotemperature energy of the crystal as a function of X. In the case X = V、differentiation of Equation (4) with res

32、pect to volume gives the equation of state:dV譏 1匸11一訐+坯窗沏qv）尹苻where*q,v) =V 9oXg,v) w(q,v) dVare the so-called Griineisen mode parameters In a perfectly harmonic crystal, phonon frequencies do not depend on the interatomic distances, hence on volume. In such a harmonic crystal Equation (5) implies t

33、hat the temperature derivative of pressure at fixed volume vanish: (dP/dT)v = 0. It follows that the thermal expansivity, p = Vl(dV/dT)P. which is given by the thermodynamical relation:(7)(9)(dp/譏(ap/av)r烏Ml=工方 3(q，v)Y(q，v"(q,v)“T qvwhere BT =V(dP/dV)T is the crystal bulk modulus, would also va

34、nish for perfectly harmonic crystals. Inspired by Equation (2), let us define Cv(q, v) = /co(q, v)nz(q, v) / V as the contribution of the v-th normal mode at the q point of the BZ to the total specific heat, and 7 as the weighted average of the various Griineisen parameters:(10)qv*q，v)G(q,v)In terms

35、 of 7, the thermal expansivity simply reads:(11)=TBrP2The vanishing of the thermal expansivity in the HA would also imply the equality of the constant-pressure and constant-volume specific heats. By imposing that the total differentials of the entropy as a function of pressure and temperature or of

36、volume and temperature coincide, and by using the Maxwell identities, one can in fact show that (Wallace 1998):(12)(13)We conclude this brief introduction to the QHA by noticing that the ansatz given by Equation (4) for the crystal free energy in terms of its (volume-dependent) vibrational frequenci

37、es gives immediate access to all the equilibrium thermal properties of the system Whether this implicit account of anharmonic effects through the volume dependence of the vibrational frequency only is sufficient to describe the relevant thermal effects, or else an explicit account of the various pho

38、non-phonon interactions is in order, instead, is a question that can only be settled by extensive numeric experienceAB INITIO PHONONSLattice dynamics from electronic-structure theorvSeveral simplified approaches exist that allow to calculate full (harmonic) phonon dispersions cc(q,v) from semi-empir

39、ical force fields or inter-atomic potentials (Briiesch 1982;Singh 1982). The accuracy of such semi-empirical models is however often limited to the physical conditions (pressure, atomic coordination, crystal structure, etc.) at which the interatomic potentials are fitted Really predictive cakulation

40、s, not biased by the experimental information used to describe inter-atomic interactions require a proper quantum-mechanical description of the chemical bonds that held matter togethe匚 This can be achieved in the framework of electronic-structure theory (Martin 2004), starting from the adiabatic or

41、Bom and Oppenheimer (BO) approximation, and using modem concepts from DFT (Hohenberg and Kohn 1964; Kohn and Sham 1965) and perturbation theory (Baroni et al. 2001).Within the BO approximation, the lattice-dynamical properties of a system are determined by the eigenvalues E and eigenfunctions of the

42、 Schrodinger equation:一弓需需+ %(R)卜(R) = E (R)Q4)where Rz is the coordinate of the /-th nucleus, Mt its mass, R indicates the set of all the nuclear coordinates, and EBO is the ground-state energy of a system of interacting electrons moving in the field of fixed nuclei, whose Hamiltonian一which acts on

43、to the electronic variables and depends parametrically upon R)reads:計 2-22|弘(R)=f刁礦話內(nèi)+糾山)+MR)(15)-e being the electron charge, VRJ(r) = -EXZ2/|r-Rz|) is the electron-nucleus interaction, and EMR)=(以/2)E冷/(Z/Zj/IRz-R)the inter-nuclear interaction energy. The equilibrium geometry of the system is dete

44、rmined by the condition that the forces acting on individual nuclei vanish:(16)whereas the vibrational frequencies, g are determined by the eigenvalues of the Hessian of the BO energy, scaled by the nuclear masses:det1 血(R)2- 3(17)Baroni etal uidd 512/9/2009 2 15 59 PMThermal Properties from Ab Init

45、io Quasi-Harnionic Phonons#Baroni etal uidd 512/9/2009 2 15 59 PMThermal Properties from Ab Initio Quasi-Harnionic Phonons#The calculation of the equilibrium geometry and vibrational properties of a system thus amounts to computing the first and second derivatives of its BO energy surface. The basic

46、 tool to accomplish this goal is the Hellmann-Feynman (HF) theorem (Hellmann 1937; Feynman 1939), which leads to the following expression for the forces:(18)Baroni etal uidd 512/9/2009 2 15 59 PMThermal Properties from Ab Initio Quasi-Harnionic Phonons#Baroni etal uidd 512/9/2009 2 15 59 PMThermal P

47、roperties from Ab Initio Quasi-Harnionic Phonons#where nR(r) is the ground-state electron charge density corresponding to the nuclear configuration R The Hessian of the BO energy surface appearing in Equation (17) is obtained by differentiating the HF forces with respect to nuclear coordinates:(1刃02

48、% (R)二dRzdRy dRy=J響警血+皿）鵲dr +鶻(20)Baroni etal uidd 512/9/2009 2 15 59 PMThermal Properties from Ab Initio Quasi-Harmonic Phonons7Equation (20) states that the calculation of the Hessian of the BO energy surfaces requires the calculation of the ground-state electron charge density, nRJ(r), as well as

49、 of its linear response to a distortion of the nuclear geometry. dnR(r)/dRj. This fundamental result was first stated in the late sixties by De Cicco and Johnson (1969) and by Pick, Cohen, and Martin (1970). The Hessian matrix is usually called the matrix of the inter-atomic farce constants (IFC). F

50、or a crystal, we can write:C；(R_R) =%(R)aa(R)a(Rz)where u5°(R) is the a-th Cartesian components of the displacement of the 5-th atom of the crystal unit cell located at lattice site R. and translational invariance shows manifestly in the dependence of the IFC matrix on R and Rz through their di

51、fference only.Density-functional perturbation theoryWe have seen that the electron-density linear response of a system determines the matrix of its IFCs, Equation (20). Let us see now how this response can be obtained from DFT. The procedure described in the following is usually referred to as densi

52、ty-functional perturbation theory (Baroni et al. 1987, 2001).(22)(23)(24)In order to simplify the notation and make the argument more general, we assume that the external potential acting on the electrons is a differentiable function of a set of parameters,九三 A/ (At= R/in the case of lattice dynamic

53、s). According to the HF theorem, the first and second derivatives of the ground-state energy read:In DFT the electron charge-density distribution, is given by:N/2nx(r) = 2Xlvi(r)|2where N is the number of electrons in the system (double degeneracy with respect to spin degrees of freedom is assumed),

54、 the single-particle orbitals, 謝入(r), satisfy the Kohn-Sham (KS) Schrodinger equation:(方2護(hù)(25)(26)H(r)=e"v"(r)and the self-consistent potential, Vsc. is given by:略=/ + F J 理纟+ BcH l(r)r-rwhere yixc is the so-called exchange-correlation (XC) potential (Kohn and Sham 1965). The electrodensit

55、y response,加入(r)/0入，appearing in Equation (23) can be evaluated by linearizing Equations (24), (25), and (26) with respect to wave-function, density, and potential variations, respectively. Linearization of Equation (24) leads to:NI2n(r) = 4%工屮；(<)屮;(D(27)M=1where the prime symbol (as in n")

56、 indicates differentiation with respect to one of the Xs. The super-script 九 has been omitted in Equation (28), as well as in any subsequent formulas where such an omission does not give rise to ambiguities. Since the external potential (both unperturbed and perturbed) is real. KS eigenfunctions can

57、 be chosen to be real, and the sign of complex conjugation, as well as the prescription to keep only the real part can be dropped in Equation (27).The variation of the KS orbitals, <|/n(r), is obtained by standard first-order perturbation theory (Messiah 1962):(H；cf= -匕)I 忙(28)where HCF = -(/i2 /

58、 2/n)(d /9r2) + Vf(r) is the unperturbed KS Hamiltonian,WscF(r)= z(r) + J K(r, r)nXr)dr(29)is the first-order correction to the self-consistent potential, Equation (26), K(r,rz) = (e2/-/!) + 6j.ixc(r)/Sn(r,) is the Hartree-plus-XC kernel, and VScfis the first ordervariation of the KS eigenvalue, zn. Equations (28-30) form a set of sclf-consistcnt equations for the perturbed system completely analogous to the KS equations in the unperturbed case Equations (24), (25), and (26)with the KS eigenvalue equation. Equation (25), being replaced by a