Ab initio up to the melting point: Anharmonicity and vacancies in aluminum
Blazej Grabowski, Lars Ismer, Tilmann Hickel, and Jörg Neugebauer
Motivation

- Figure 1: High-temperature isobaric heat capacity of aluminum from experiment and ab initio calculations. See Ref. [5] for the original references of the experimental data. The melting temperature of aluminum, 933 K, is denoted by Tm and the Boltzmann constant by kB.
At elevated temperatures, the heat capacity of metals significantly deviates from the behavior predicted within the perfect harmonic lattice model. This was pointed out already in 1921 in a seminal work by Born and Brody [1]. In this and in many subsequent studies, various mechanisms have been proposed to explain these deviations. Only recently, our ab initio calculations for a wide range of non-magnetic metals showed that the predominant part of these deviations can be explained on the basis of quasiharmonic (going beyond linear Grüneisen theory) and electronic excitations (see Ref. [2] and the discussion here). A generalization to magnetic systems was given in Ref. [3] and is also discussed here. However, the subtle balance between further contributions, in particular explicit anharmonicity (see below for definition) and vacancies, could not be clarified yet even for simple elementary metals. [4]
A major challenge in determining the balance between the explicit anharmonicity and vacancy contribution to the high-temperature heat capacity has been up to now a significant scatter in experimental data. A striking example is the heat capacity of aluminum shown in Fig. 1. The ambiguous experimental situation has considerably hampered empirical approaches in determining the contributions and led to a long standing debate whether explicit anharmonicity or vacancies dominate the heat capacity of aluminum. (See Ref. [5] for a detailed discussion.)
Ab initio methods, i.e., methods based solely on quantum mechanical principles, are well suited to tackle such a challenge. However, a direct ab initio calculation of the explicit anharmonicity contribution to the heat capacity has not been possible up to now due to large computational requirements. In the present study [5], we have developed and applied numerically and computationally highly efficient ab initio methods based on density-functional theory to sample the configuration space. With these methods we were able to calculate the heat capacity of aluminum including all relevant excitation mechanisms. We particularly focused on the reduction of controllable errors to less than 1 meV/atom in the free energy.
Theory
Any thermodynamic property such as the heat capacity can be obtained from the free energy surface F(V,T), where V is the volume and T the temperature. Since aluminum is a non-magnetic element F(V,T) can be decomposed into:
Here, Fel is the electronic free energy including the T=0 K binding energy, F qh+Fah=Fvib is the vibronic free energy Fvib due to the ionic motion consisting of the quasiharmonic free energy Fqh and explicitly anharmonic free energy Fah, and Fvac is the free energy contribution due to vacancies. This decomposition is justified by the fact that the various excitation mechanisms reside on different time scales and therefore their coupling is small compared to the individual contributions given in Eq. (1). The evaluation of Fel and Fqh is discussed in detail here. Fah is the remaining part of Fvib which is not captured by the quasiharmonic approximation. It therefore includes the explicit temperature dependence of the phonons, i.e., phonon-phonon coupling. A common approach, which we have also applied in this study, is the thermodynamic integration [6]. Fvac is usually calculated based on the assumption of a constant vacancy supercell. We found that one needs to go beyond this to achieve the desired accuracy. We have therefore derived a set of coupled equations, which allow to treat Fvac based on a self-consistently optimized vacancy supercell. [5]
Methodology

- Figure 2: Logarithmic plot of the single processor CPU time needed to calculate an ab initio anharmonic free energy (for the system described in the text) using direct molecular dynamics (MD), thermodynamic integration (TI), and the proposed method (UP-TILD = UPsampled Thermodynamic Integration using Langevin Dynamics).
The main challenge in calculating Fah from ab initio arises from long computational times. A direct statistical sampling of the configurational space using molecular dynamics (MD) would require ≈107 MD steps. Taking the quasiharmonic solution in a thermodynamic integration as a reference, this number can be reduced to ≈104. This is, however, still to large for a direct ab initio sampling, since one MD step typically requires about 1 CPU hour for the considered system which results in calculation times of several thousand CPU hours.
In this project, we developed a hierarchical scheme which allows to reduce the number of expensive ab initio calculations to ≈102 (cf. Fig. 2). The basic idea (see Ref. [5] for details) is to use density-functional theory calculations (e.g. low cutoff and k-points) to sample the configurational space and to determine Fah,low. In a next step, a special set of M MD structures is chosen and the difference ΔE between a low and high converged energy, Elow and Ehigh respectively, is calculated for each of these structures. The average of ΔE over all special structures
is then used to obtain the desired highly accurate anharmonic free energy from:
Using this scheme, we did not only calculate Fah for the perfect bulk, but we were also able to compute for the first time a fully ab initio Fvac(V,T) surface including all relevant excitations:
Results and conclusions
A selection of our ab initio results regarding the balance between the various physical mechanisms contributing to the isobaric heat capacity of aluminum is given in Fig. 3.

- Figure 3: Ab initio calculated contributions to the isobaric heat capacity at zero pressure of aluminum in units of the Boltzmann constant kB. The melting temperature Tm is indicated by the vertical dashed line. a) Contributions arising from electronic, anharmonic, and vacancy (vac) excitations. The sum of these contributions is given by the black solid line. The right axis are scaled with respect to the 'full' heat capacity from b) at Tm. b) The quasiharmonic contribution and the full heat capacity corresponding to the sum of the quasiharmonic one and the contribution given by the black line in a).
The major contribution is given by the quasiharmonic contribution (Fig. 3b), whereas the other contributions (Fig. 3 a) can be regarded as small perturbations. A crucial observation is that the anharmonic contribution is negative and in fact cancels the positive electronic contribution over a wide range of temperatures. The finding of a negative anharmonic heat capacity is in contrast to the common belief according to which it was either vanishing or positive. [7,8] The vacancies, on the other hand, yield a positive non-linear increasing contribution, which is however not as strong as the non-linear increase in the experimental isobaric heat capacity observed in several measurements. This fact is demonstrated in Fig. 1 where we compare our 'full', i.e., containing all contributions listed in Eq. (1), ab initio isobaric heat capacity to an extensive set of experiments. For further information regarding the obtained results and developed methods the reader is referred to Ref. [5]. Some of the issues presented in Ref. [5] are: results for two exchange-correlation functionals (LDA and GGA), results for further thermodynamic quantites such as vacancy concentration and expansion coefficient, discussion of the influence of the various free energy contributions on the vacancy quantities, and development of an intuitive approximate model to describe the anharmonic free energy.
We believe that the achievable numerical accuracy of our approach will open a number of interesting research topics. This accuracy is crucial to systematically check the performance of density-functional theory with respect to its ability in describing finite temperature material properties.
References
[1] M. Born and E. Brody, Z. Phys. 6, 132 (1921).
[2] B. Grabowski, T. Hickel, and J. Neugebauer, Phys. Rev. B 76, 024309 (2007).
[3] F. Körmann, A. Dick, B. Grabowski, B. Hallstedt, T. Hickel, and J. Neugebauer, Phys. Rev. B 78, 033102 (2008).
[4] Y. Kraftmakher, Phys. Rep. 299, 79 (1998).
[5] B. Grabowski, L. Ismer, T. Hickel, and J. Neugebauer, submitted to Phys. Rev. B.
[6] L. Ismer, Ph.D. thesis, University of Paderborn (2008).
[7] C. R. Brooks and R. E. Bingham, J. Phys. Chem. Solids 29, 1553 (1968).
[8] R. C. Shukla, C. A. Plint, and D. A. Ditmars, Int. J. Thermophys. 6, 517 (1985).