Development of a tight-binding potential for bcc-Zr. Application to the study of vibrational properties

We present a tight-binding potential based on the moment expansion of the density of states, which includes up to the fifth moment. The potential is fitted to bcc and hcp Zr and it is applied to the computation of vibrational properties of bcc-Zr. In particular, we compute the isothermal elastic constants in the temperature range 1200K<T<2000K by means of standard Monte Carlo simulation techniques. The agreement with experimental results is satisfactory, especially in the case of the stability of the lattice with respect to the shear associated with C'. However, the temperature decrease of the Cauchy pressure is not reproduced. The T=0K phonon frequencies of bcc-Zr are also computed. The potential predicts several instabilities of the bcc structure, and a crossing of the longitudinal and transverse modes in the (001) direction. This is in agreement with recent ab initio calculations in Sc, Ti, Hf, and La.


I. INTRODUCTION
A fundamental problem in Condensed Matter Physics is how to project microscopic interactions in many-body systems onto a description in terms of a reasonably small number of degrees of freedom [1]. This is particularly important in computer simulation studies. Atomic-level simulations are, intrinsically, large-scale calculations and in spite of the huge improvement nowadays available in computing capacity, an efficient method to rapidly evaluate energies which also treat forces in a physically realistic way, still remains a major difficulty to be solved.
The computation of solid properties always requires a particular choice for the interatomic potential. In the choice there is a compromise between physics and efficiency. Computational efficiency makes empirical or semiempirical potentials desirable, but at the same time one should require that the underlying physics behind the model potential is able to reproduce the properties of the system or at least those of interest.
Several years ago Friedel [2] suggested that the starting point in the description of transition metals (TM) is a band-picture with a strong d-character. In this sense, the remarkable parabola-like behaviour of the cohesive energy and the bulk modulus exhibited by most TM as a function of the number of d-electrons [2,3] clearly indicates that cohesion is mainly dominated by the d-states. This has motivated the development of many-body potentials based on the tight-binding (TB) approximation [4,5] (For a review see Ref. [6]). They are semiempirical in nature, which makes them very appealing for computer simulation studies and at the same time they incorporate the band character of the metallic cohesion so that the attractive part of the interatomic energy turns out to be many body.
In the present work, we develop a TB potential based on the moment expansion of the density of states, which includes up to the fifth moment. It has been suggested that this is the lowest order needed to reproduce the general trends in the elastic constants of hcp and fcc TM as a function of band filling [7].
Zirconium, as is the case of other metals and alloys, is close packed at low temperatures, but undergoes a structural phase transition of the Martensitic type to a bcc structure at higher temperatures [8,9,10]. Some features of this phase transformation in Zr were previously studied by using the second moment approximation TB model [11,12]. Here, we shall not focus on the structural phase transition itself but rather concentrate our interest on the elastic properties of the high-temperature bcc-phase, which is accepted to be mainly stabilized by entropy effects [13,14,15]. In particular we calculate the temperature behaviour of the relevant elastic constants by using standard Monte Carlo simulation techniques. The results thus obtained are compared with the available experimental data. The agreement is rather satisfactory, but it does not allow us to draw conclusions about the reliability of the potential. This is provided by the computation of the phonon dispersion curves for the bcc-Zr at T = 0K. Indeed, the comparison with recent ab initio calculations [16] is indicative that most of the fundamental physics governing the vibrational properties of bcc-Zr is contained in the interatomic potential model. The paper is organized as follows. The next section is devoted to the construction of the interatomic potential. In Sec. III we describe the fitting procedure to Zr. In Sec. IV we present and discuss the results and in Sec. V conclusions are drawn.

II. DEVELOPMENT OF THE INTERATOMIC POTENTIAL
The interatomic potential is developed following the tight binding bond model (TBBM) by A. P. Sutton et al. [17,18] in the two-centre orthogonal approximation. The basis set of the TB Hamiltonian only includes d atomic orbitals and the crystal field interactions are neglected.
The cohesive energy is decomposed into two terms, The term E pair is an empirical pair potential which stands for electrostatic and exchange-correlation interactions [17], and the bond energy is where E F is the Fermi energy, ǫ i are the on-site Hamiltonian matrix elements in the atomic orbital representation, and n i (E) are the local densities of states (LDOS). In this TB model, the s-d hybridization, which is known to make a considerable contribution to the cohesive energy of TM [19], is not treated appropriately. Nevertheless, it can be assumed that such contribution is included implicitly in the pair term, as in the paper by Girshick et al. [20], or either that this contribution should be proportional to the d band width and therefore it is included in the bond term.
The condition of local charge neutrality is fulfilled by defining local Fermi energies. This is done through the relation, where N d is the number of d electrons per site and the atomic orbital energy is chosen to be the energy zero, ǫ i = 0 ∀i. This method is equivalent to shifting the LDOS rigidly, whereas adjusting the on-site energies selfconsistently, as proposed by Sutton et al., means that the shifts of the LDOS are accompanied by a distortion of their shape. The LDOS are constructed from their moments, µ (n) i , following the formalism of the recursion method of Haydock [21] as if the moments corresponded to an s band. In this way the computed LDOS are rotationally invariant [22]. That is, from the second to the fifth moments we compute the coefficients b 1 , a 1 , b 2 and a 2 of the recursion method. Since the coefficients a n , b n are convergent oscillating series [23] their limit a ∞ , b ∞ is estimated as a ∞ = (a 1 + a 2 )/2 and b ∞ = (b 1 + b 2 )/2, and the coefficients a n , b n with n > 2 are assumed to be equal to a ∞ , b ∞ which gives rise to the well known square root terminator of the continued fraction expansion of the diagonal elements of the Green function [24]. The integration of the LDOS in order to obtain the bond energy is performed numerically.
To completely specify the bond energy we need to define the functional form of the off-diagonal Hamiltonian matrix elements in the atomic orbital representation, known as coupling strengths or hopping integrals. In the two-centre approximation the coupling strengths between sites i and j can be written as where Φ αβ (r ij /r ij ) is the angular dependence given by the symmetry of the d orbitals, as shown by Slater and Koster [25], and R(r ij ) is the dependence on distance of the ddσ, ddπ and ddδ bonds, which is assumed to be equal for all of them. The relative strength of the couplings V ddσ : V ddπ : V ddδ are the canonical values -6 : 4 : -1.
The model has, in principle, six fitting parameters: A, ξ, p 1 , p 2 , V ddδ and q. Nevertheless, the number of electrons, N d , and the cutoff distances r 1 , r 2 ,r 1 , andr 2 are also used as fitting parameters although they are not completely free since they are constrained by their physical character.

III. FITTING TO ZR
In this section we discuss the procedure used to fit the model parameters to the zero temperature properties of Zr, both in the hcp and the unstable bcc phases.
Since the bcc structure is only stable at T > 1135K, the values of some properties used in the fitting procedure are extrapolated from the high temperature experimental data, whereas others are obtained from ab initio calculations at T = 0K. Most of the ab initio results have been calculated in this work using the WIEN97 code [26], which follows the Full Potential Linearized Augmented Plane Wave method within density functional theory. All the computed quantities are obtained in the generalized gradient approximation (GGA) by Perdew, Burke and Ernzerhof [27] for the exchange-correlation potential.
The properties of bcc-Zr that we have considered in the fitting procedure are: the cohesive energy, the lattice parameter, the elastic constants and the unrelaxed vacancy formation energy. Moreover, we have also taken some properties of hcp-Zr at T = 0K into account, mainly its elastic constants, lattice parameter, the c/a ratio, and the energy difference between both structures.

A. Fitting of the bond energy
The coupling strengths are essentially short-range functions. Therefore, we impose that they fall to zero between the second and third nearest neighbours and choose r 2 = 1.24a bcc . For the matching point of the tail of the function we chose r 1 = 0.92a bcc , which lies between the first and second nearest neighbours. Notice that in this way the tail of the coupling strengths interferes with the fitting procedure. The values chosen for these parameters are those which allow a better reproduction of the properties of hcp-Zr.
An important property of the elastic constants exploited in the fitting procedure is that the Cauchy pressure given by an interatomic potential which satisfies the mechanical equilibrium conditions only depends on the many-body term of the potential, that is, the Cauchy pressure is independent of the pair term [20]. Considering both the bcc and the hcp structures we have three Cauchy pressures to deal with. One of these, the Cauchy pressure C 12 − C 66 of the hcp lattice is affected by an internal relaxation of the lattice under strain. In order to avoid the determination of the internal relaxation and its effect on the Cauchy pressure during the fitting procedure, as the test value for this quantity, we use only its homogeneous contribution (without internal relaxation) obtained from the experimental Cauchy pressure by subtracting the inhomogeneous contribution determined from ab initio calculations. In the hcp lattice the elastic constants affected by internal relaxations are C 11 , C 12 and C 66 = (C 11 − C 12 )/2, and the inhomogeneous contribution to C 12 is the same as that of C 11 but with the opposite sign. Therefore, we need only to estimate the inhomogeneous contribution to C 11 which is I ≃ 7 GPa [28]. Using the experimental values of the elastic constants together with the ab initio results, the three Cauchy pressures are, (3.1) Another property which strongly depends on the many-body term is the unrelaxed vacancy formation energy, (E f V ) unrelx , since the contribution of the pair potential term to this quantity is equal to its contribution to the cohesive energy. For the present interatomic potential, in the bcc lattice we get the following approximate relation, Taking into account (2.1) we therefore get, This result is consistent with the value predicted by the renormalized-atom method of Gelatt et al. [19] provided that the contribution to the cohesive energy of the s-d hybridisation is assumed to be contained in the bond term From the literature we obtain physical values for the number of electrons, N d = 2.536 [29], band width, W =7.8 eV [30], and the derivative of the band width with respect to the atomic volume, −3Ωd ln W/dΩ = 3.97 [30,31], which is related to the parameter q. In order to extract the value of the parameter q from this relation we approximate the band width as W ∝ µ (2) , that is, proportional to the square root of the second moment of the density of states [3], and using the radial form of the coupling strengths defined in Eq. (2.5), we obtain q ≃ 4.36a −1 bcc . In the present model the band width is W = E top − E bottom = 4b ∞ which is proportional to the parameter V ddδ . Therefore we can use W instead of V ddδ as the fitting parameter without loss of generality.
Using these physical values for the parametres N d , W and q, in the bcc lattice the interatomic potential gives a value of E bond ≃ −7.5 eV, which is consistent with the value predicted by the unrelaxed vacancy formation energy (Eq. (3.3)). Nevertheless, the Cauchy pressure obtained with these parameters, (C 12 − C 44 ) bcc ≃ 104 GPa, is unacceptably high (see Eq. (3.1)). We have, therefore, to find an alternative way of fitting them. We shall proceed by paying more attention to the directly observable physical quantities such as cohesive energy, lattice parameter, elastic constants and vacancy formation energy rather than to the physical interpretation of the parameters.
For a given value of N d , the parameters W and q are determined numerically to reproduce the Cauchy pressure of the bcc lattice given in Eq. (3.1) and the bond energy, E bond = −8.0eV (Eq. (3.3)). We repeat this fitting procedure for different values of N d and compute the Cauchy pressures of the hcp lattice. The results are given in Table I. These Cauchy pressures are computed using the lattice parameter and c/a ratio of the hcp lattice at which this structure is stabilized by the interatomic potential. For this purpose we use an arbitrary pair contribution fitted to the lattice parameter, cohesive energy and elastic constants of the bcc structure. In all cases, a hcp and c/a are close to the experimental and ideal values respectively, although we observe a tendency of c/a to increase as the number of electrons is increased. The value of N d which allows reproduction of the three Cauchy pressures simultaneously (Eq. (3.1)) is N d = 1.45, which is substantially lower than its physical value, N d = 2.536 [29]. Nevertheless, there are two additional reasons for adopting the value of N d = 1.45, i) TB studies on the relative stability of hcp and fcc lattices as a function of N d [29] have shown that the hcp lattice is only stable for N d < 2 (in the range 0.5 < N d < 4.5), and ii) the low value of the c/a ratio in group IV hcp metals can be justified in terms of the different contribution to the bond energy from nearest neighbours located on the basal plane and those located off it only if N d < 2 [17,32].
The radial dependence of the coupling strengths is shown in Fig. 1.

B. Fitting of the pair potential
The cutoff distance of the pair potential is chosen to ber 2 = 2.211a bcc which lies between the 7th and 8th nearest neighbours in the bcc structure. This value has been adjusted in order to reproduce the properties of the hcp lattice. The value ofr 1 is almost irrelevant since it does not really affect the functional form of the pair potential. We chosẽ r 1 = 1.2a bcc . For given values of the parameters p 1 , p 2 and ξ, the parameter A is determined from the mechanical equilibrium condition of the bcc lattice The parameters p 1 , p 2 and ξ are given by the cohesive energy, and elastic constants C 11 and C 12 of the bcc lattice. Notice that the elastic constant C 44 is automatically reproduced by fitting C 12 since the Cauchy pressure is fixed by the bond energy.
The functional form of the pair potential is shown in Fig. 1.
The values of the fitted parameters are given in Table II and the properties of bcc and hcp Zr used for the fit are given in Tables III and IV respectively. In the fitting procedure we have priorized the reproduction of the properties of bcc-Zr since the interatomic potential will be used to study the vibrational properties of this structure. This have led to different accuracy in the fitted quantities of bcc and hcp Zr (see Tables III and IV).

IV. CALCULATION OF THE VIBRATIONAL PROPERTIES OF BCC-ZR
In this section we use the interatomic potential already obtained to compute the phonon dispersion curves of the bcc lattice at T = 0K, and the elastic constants in the temperature range 1200K < T < 2000K, where experimentally, the bcc structure is stable.
A. Phonon frequencies of bcc-Zr at zero temperature In order to obtain the phonon frequencies of bcc-Zr, we first compute the interatomic force constants by means of standard numerical differentiation. It is worth pointing out that the interatomic potential gives rise to long-range force constants. In particular, the range of the pair potential contribution is the range of the potential, which in the present case is up to the 7th nearest neighbours. The range of the bond term contribution is much larger. For coupling strengths extending up to the second nearest neighbours and including up to the fifth moment of the LDOS in the computation of the bond energy, the range of the force constants extends up to the 22th nearest neighbours. This is due to the many-body character of the bond energy together with the dependence of the high order moments on the position of distant atoms.
In Table V we show the results for the computed force constants. From these values we construct the dynamical matrix and compute the phonon frequencies in the harmonic approximation along the high symmetry directions of the reciprocal space (Fig. 2). Since the interatomic potential is fitted to the T = 0K elastic constants of bcc-Zr, the slope of the phonon dispersion curves around the Γ point is expected to be correct. There are several features of the phonon dispersion curves that deserve special comment: 1. The whole T1(ξ ξ 0) branch is unstable.
2. The T1(ξ ξ 2ξ) branch has a positive slope around the Γ point (consistent with the associated combination of elastic constants), but it rapidly softens and becomes unstable. Before matching the T1(ξ ξ 0) branch at the N point, it becomes stable again at about ξ = 1/3.
3. The softening of the L(ξ ξ ξ) branch around ξ = 2/3 observed experimentally at high temperature is, at zero temperature, an instability which gives rise to the ω phase. The frequency of the ξ = 2/3 phonon is approximately zero, as was obtained by K.M. Ho et al. [33] from ab initio calculations, and the minimum of the branch is at about ξ = 17/24 4. In the (001) direction there is a crossing between the longitudinal and transverse branches.
Several of these features, not observed experimentally at high temperature, were obtained from ab initio calculations in Sc, Ti, Hf, and La at T = 0K [16]. In these materials the whole T1(ξ ξ 0) branch has imaginary frequencies (the phonons in the (112) direction have not been computed). The instability around the ξ = 2/3 L(ξ ξ ξ) mode is also predicted, although the minimum of the branch is not located at ξ = 17/24 but at ξ = 7/12 for all elements (except Sc, which has the minimum at ξ = 2/3). Finally, the crossing of the longitudinal and transverse branches in the (100) direction is also observed.
Moreover, in these materials the T(ξ ξ ξ) branch has imaginary frequencies around the Γ point. The slope of this branch is given by the elastic constant C 11 − C 12 + C 44 , which in Zr is positive, and thus this feature cannot be observed.
From the force constants we also compute the elastic constants following the method of long waves [34], and recover the values obtained by means of homogeneous deformations. This can be used as a test of the internal consistency of the interatomic potential [17].

B. High temperature elastic constants of bcc-Zr
The bcc-Zr high-temperature elastic constants are obtained from Monte Carlo (MC) simulations in the canonical ensemble, (T, V, N ), using the atomic volume obtained from MC simulations in the isobaric-isothermal ensemble, (T, P, N ), at zero pressure. The theoretical background of such simulations can be found in the work by I. R.
The second-order isothermal elastic constants are computed using the fluctuation formula [36] where F is the Helmholtz free energy, ǫ ij are the elastic strains, V is the total volume, k B is the Boltzmann constant, T is the temperature, and N is the number of atoms. The derivatives of the cohesive energy with respect to the elastic strains are computed numerically. Since the second derivative of both the pair potential and the coupling strengths is not continuous at the cutoff distance, the numerical derivative must be calculated using the same radial and pair functions in both the strained and the unstrained state of the lattice, for each of the atomic pairs. That is, if for a given pair of atoms the distance, r, is r 1 < r < r 2 , the radial function used for the computation of the bond energy of both the strained and unstrained lattice will be that defined in this region, regardless of the fact that in the strained lattice we may have r > r 2 .
The simulations are performed on a 4 × 4 × 4 × 2 = 128-site bcc lattice with periodic boundary conditions. The attempted configurational changes are single atom displacements and after each N of these attempts (= 1 Monte Carlo step (MCS)), in the isobaric-isothermal ensemble a volume change is also proposed. The simulations are 10 5 and 5 · 10 4 MCS long in the isobaric-isothermal and canonical ensemble respectively.
Due to the many-body character of the bond energy, the change in the total cohesive energy due to a singleatom movement involves the recalculation of the contribution to bond energy of about 65-110 atoms, depending on temperature. Nevertheless, this computation can be highly optimized and only requires about six times the CPU time needed to compute the contribution to the bond energy of a single atom.
In Fig. 3 we show the computed lattice parameter vs. temperature. Although at temperatures above 1200K the computed lattice parameter is about 0.6% smaller than the experimental value, the thermal expansion coefficient (slope) is reproduced to great accuracy β = 1/V (∂V /∂T ) = 3.0 · 10 −5 K −1 (β exp = 2.8 · 10 −5 K −1 [37]). Notice that the linear extrapolation of the computed lattice parameter to zero temperature does not match the fitted value. This is because at low temperature the thermal expansion given by the interatomic potential is strongly nonlinear. This is probably related to the fact that the bcc structure is unstable, although this behavior is not observed when comparing the high temperature experimental results [37] with the zero temperature ab initio calculations (see Table III).
In Fig. 4 we show the temperature dependence of the relevant elastic constants obtained from Monte Carlo simulations. The main success of the present interatomic potential is that it renders a positive C ′ at high temperature. Moreover, the value predicted for the whole set of elastic constants, C 11 , C 12 , and C 44 is rather accurate. The main failure is that it is unable to reproduce the large value of C ′ observed experimentally. This failure comes mainly from the values of the elastic constant C 12 which experiments have shown to decrease strongly with temperature. This marked decrease of C 12 , together with the nearly constant behaviour of C 44 , means that the Cauchy pressure decreases with temperature. We were unable to reproduce such behaviour. In several tests during the parametrization of the interatomic potential we always found a Cauchy pressure nearly independent on temperature.
In Table VI we show separately the different contributions to the elastic constants: the Born term, the fluctuation term, and the kinetic term. We get the rather unusual result that the fluctuation term of C 12 is nearly zero. This means that the strains in different directions are uncorrelated in the simulations Since the contribution of the fluctuation term to C 12 is usually negative and the interatomic potential is unable to reproduce the low value of C 12 observed experimentally, we conclude that this lack of correlation given by the interatomic potential is possibly unphysical.

V. DISCUSSION AND CONCLUSIONS
In the present paper we have developed a TB interatomic potential suitable for the study of the vibrational properties of bcc Zr. The interatomic potential has been fitted to the T = 0K properties of Zr in both the hcp and the bcc structures. Although among the vibrational properties, only the elastic constants are used in the fitting procedure, the TB potential shows a remarkable capacity of predicting the T = 0K phonon frequencies of the bcc structure along the high symmetry directions studied. As regards the high-temperature elastic constants, the general trends are reproduced, especially the stability of the bcc structure with respect to the shear associated with the elastic constant C ′ . Nevertheless, the interatomic potential is unable to reproduce the temperature decrease of the Cauchy pressure.
The reliability of the experimental values of the high-temperature elastic constants should, however, be questioned. The elastic constants cannot be obtained from measurements of the velocity of acoustic waves in the material because the temperature at which the bcc phase is stable is too high. Heiming et al. [37] therefore derived the elastic constants from the force constants obtained from a fit to the phonon dispersion relations. In order to keep the number of force constants reasonably small, in the fitting procedure they impose that the range is up to the fifth neighbour shell, which is rather short. On the other hand, the elastic constants obtained depend critically on the frequencies of the phonons close to the Γ point of the Brioullin zone, and thus have large error bars. We have tried to derive the elastic constants from the phonon frequencies of Heiming et al. [37], but we only were able to reproduce the elastic constant C ′ to any accuracy. The values of all the other elastic constants strongly depend on the phonons considered and the method of fitting.
In spite of the remarkable success of the interatomic potential in reproducing the T = 0K phonon frequencies, we should mention the difficulties we have found during the fitting procedure, and discuss which features of the TB potential are expected to correctly reproduce the physics of the material and which are not.
The first point concerns the range of the hopping integrals, which in fact is too small. In the hcp lattice they should fall between the second and third nearest neighbours, at least, and only the nearest neighbours are taken into account. This leads to a bond energy contribution in the hcp lattice that is smaller than in the bcc lattice, which is just the opposite to the expected result. This fact is compensated by the pair potential contribution to give the correct energy difference between both structures, but it is still clearly reflected in the unrelaxed vacancy formation energies.
The reason for choosing such a low value for the cutoff distance of the hopping integrals is computational convenience. The range of the hopping integrals in the bcc lattice is up to second nearest neighbours (14 atoms involved). This means that in the perfect lattice at T = 0K the number of neighbours involved in the computation of the moments is 64 (up to the sixth coordination shell). Nevertheless, at T = 2000K the atoms are far from their equilibrium positions and this rapidly increases the number of neighbours involved in the computation of the moments up to ≃ 110. In order to correctly compute the moments we must therefore take into account the neighbours up to the 13th coordination shell (258 atoms). An increase in the cutoff distance of the hopping integrals involves an increase of the number of neighbours to be taken into account in the computation of the moments, and thus, we decided to choose the lowest possible value which allowed us to obtain physically reasonable results.
The second point concerns the pair potential contribution. During the fitting procedure we observed that the capacity of the interatomic potential to simultaneously reproduce the properties of both the hcp and the bcc structure is strongly dependent on the range of the pair potential. In other words, if we take a different range to that used in this paper, the results are rather worse. This is indicative that the geometry of the different coordination shells has an important effect on the elastic constants. Moreover, although at the cutoff distance the pair potential and its first derivative are continuous, the decay to zero is still sharp, and the contribution to the elastic constants by the last coordination shell is unphysically high.
Finally, we should mention that the s-d hybridization, which is not explicitly included in the TB potential, has an important contribution to the cohesive energy (about -2eV [19]). We have considered only d atomic orbitals in the basis set in order to minimize computation time and the complexity of the TB potential.
The problems encountered when trying to use the physical values for the quantities N d , W , and q have already been discussed. Nevertheless, the treatment of these quantities as fitting parameters gives enough flexibility to the interatomic potential to reproduce to reasonable accuracy all the magnitudes described in the paper.
A significant improvement in the behaviour of the elastic constants requires a better determination of the Fermi energy, together with a more detailed description of the DOS, especially around this point. Nevertheless, the inclusion of high order moments into the interatomic potential is computationally very expensive.