A Generalized Approach Obeying the Third Law of Thermodynamics for the Expression of Lattice Stability and Compound Energy: A Case Study of Unary Aluminum

Recently, Hillert and Selleby proposed a simple method for expression of the lattice stability or Gibbs energy of formation that does not violate the third law of thermodynamics. This method describes the derivation of the Gibbs energy function from high temperatures down to 0 K by interpolation, instead of extrapolation from room temperature to 0 K. In the present work, their original method is discussed in terms of determination of the characteristic parameter values. Keeping the essential interpolation character of their method, a generalized approach is presented for expressing the lattice stability through parameter optimizations. This approach retains the zero point entropy of substances and is in line with the development of the third generation CALPHAD databases. Using the Al unary system as a case study, the lattice stabilities of the hcp and bcc phases are investigated. The respective Einstein temperatures are also evaluated. At high temperatures, the present descriptions reproduce the lattice stabilities suggested by SGTE for the existing second generation of databases, with a reasonable accuracy. More importantly, information from ab initio calculations (total energy at 0 K) is also used for this optimization and the present method results in a physically sounder description of thermodynamic properties at lower temperatures down to 0 K. The present approach provides a simple and flexible way to estimate the lattice stabilities, with potential applicability for the Gibbs energy of formation of stoichiometric compounds and the excess energy of solution phases, in accordance with the third law of thermodynamics.


Introduction
In the CALculation of PHAse Diagrams (CALPHAD) method, the Gibbs energy is described as a function of temperature, pressure and content of components. The parameters of the Gibbs energy function based on some model with a mathematical expression are determined so as to represent the information of thermodynamic quantities and phase diagram. The description of the pure elements by Scientific Group Thermodata Europe (SGTE) using the empirical polynomial functions [1] has widely been accepted, although it is restricted at temperatures mainly down to 298.15 K and the main application is in the high-temperature region. Until now, a wide range of thermodynamic databases has been developed, based on these unary descriptions, and are known as the second generation of thermodynamic databases.
Conversely, the third generation of CALPHAD databases with more physical basis is being developed, aiming for the thermodynamic descriptions reliable from 0 K up to high temperatures above the melting temperature. The ideas of the thermodynamic models considering each physical effect, including the Debye or Einstein model, were proposed at the Ringberg workshop, [2][3][4] and later some attempts were made to develop an appropriate model. [5][6][7][8] In the case of unary Fe studied by Chen and Sundman, [5] the specific heat was described with help from the model shown later as Eq 7, i.e., the harmonic lattice vibration described by the Einstein model, electronic excitation, anharmonic lattice vibration and magnetic contribution were modeled separately. More efforts have been made to assess the thermodynamic parameters based on this model using experimental information including specific heat and enthalpy in unaries Fe, [5,9] Mn [10] and Co. [11] One of the challenges in developing the thermodynamic database is to reasonably estimate the lattice stabilities of metastable and unstable phases, i.e., the differences in Gibbs energy between stable and metastable (unstable) phases for an element. The pressure dependence of Gibbs energy is not discussed in this paper. The temperature dependence of the Gibbs energy of a metastable (unstable) b phase, relative to that of a stable a phase is often described using the following approximation as a linear expression for many pure elements in the SGTE database. This simple approach has been adopted due to the lack of experimental information [1,12] : This expression indicates that both phases have the same specific heat and that the difference in entropy of the b phase against the a phase is b, independent of temperature. Therefore, the entropy of the b phase extrapolated to 0 K is not zero, which is in disagreement with the third law of thermodynamics. The situation is often similar for the Gibbs energy of formation in compounds and excess energy of solution phases. Alternatively, one simple way to extend the Gibbs energy function and their derivations to 0 K by interpolation instead of extrapolation has been recently proposed by Hillert and Selleby [13] as Eq 2: During the last decade, the methodology of the theoretical calculations and computer technique development have made the ab initio calculations increasingly applied as a powerful tool for calculating the thermodynamic properties of materials, particularly for phases that cannot be experimentally examined. For example, in the case of Al, the stable fcc phase has been studied by both experiments and theories. [14][15][16][17][18][19][20][21] It was extensively discussed in terms of the individual contributions to the heat capacity (or entropy). [7,[22][23][24] The largest contribution is, of course, the harmonic phonon which can be described by Debye or Einstein models. The anharmonic effects are also important despite its smaller amount of contributions. One of them is the phonon frequencies shift due to the thermal expansion, often described by the quasiharmonic model. This effect occupies a dominant part of heat capacity beyond the harmonic contribution. The other is referred as the explicitly (or high-order) anharmonic effect. This contribution is negative in the fcc Al, which is not a general result in substances, and partially compensates the electronic excitations. [23,24] The vacancy contribution is small. [24] The hcp and bcc phases are metastable and unstable at 0 K under ambient pressure in Al, respectively. [25][26][27][28] The stability of these phases has been calculated and the effect of pressure has been a main issue because these structures become stable at high pressures. [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] From a practical viewpoint of CALPHAD, the Gibbs energy has to be defined/estimated even for unstable phases at ambient pressure (such as bcc Al), which determines the applicability and quality of the extrapolation of the CALPHAD data into high order systems (e.g. Al solubility in other metals with bcc structure) or higher pressures (of pure Al). Moreover, reported experimental information of thermodynamic quantities at finite temperatures tends to be limited compared to that from groundstate calculations; nevertheless, it is necessary to assess the lattice stability for the application of CALPHAD data in multicomponent systems. One effective way to assess the lattice stability valid from 0 K up to high temperatures beyond liquidus is to couple the ab initio calculations with the existing SGTE descriptions into a new thermodynamic model implemented in the third generation CALPHAD database. The above mentioned Eq 2 is helpful to describe the lattice stability in agreement with the third law of thermodynamics while using both the SGTE description for T [ 298.15 K and ab initio calculations at T = 0 K.
It should be emphasized that, for the sake of zero-point vibrational entropy, Eq 2 is applicable on the temperaturedependence of not only lattice stability, but also the formation energy of stoichiometric compounds and excess energy of solution phases. The abundant well assessed thermodynamic parameters in the existing second generation CAL-PHAD are valuable. During the development of third generation one may start the assessment of a new chemical system with direct parameter transition into the format of third generation description by using Eq 2. It is worthy to discuss and demonstrate how to apply Eq 2 with a real system i.e. the unary Al in different cases where the phase is metastable or unstable. In the original paper by Hillert and Selleby [13] the choice of characteristic values T 1 and n in Eq 2 has not been sufficiently discussed. In this paper, we examine the applicability of the expression of Eq 2. First, we test various T 1 and n and then optimize the values using a -bT at high temperatures and ab initio results at 0 K. As case studies, the lattice stabilities of the hcp and bcc Al were chosen and the information of a -bT was taken from the SGTE database. From these trials, we propose a modified approach to express the lattice stability, Gibbs energy of formation of stoichiometric compounds and excess energy of solution phases from 0 K to high temperatures that were described by a linear expression in existing database.
In summary, the present paper deals with a methodology on the development of third generation CALPHAD database coupling the information at high temperatures from the second generation and 0 K from ab initio calculations. The utilized modelling quantity is Gibbs energy. Many efforts have been made to bring in more physical factors into such thermodynamic modelling in the development of third generation CALPHAD database. However, it should be clarified that there is neither any ambitions/attempts nor feasible to develop a pure physical model for Gibbs energy descriptions and consequently phase equilibria predictions. To better demonstrate the methodology meanwhile not lost the generality, the case study on unary Al was presented here for metastable hcp and unstable bcc.

The Generalized Lattice Stability
Hillert and Selleby [13] have proposed Eq 2 and also suggested the combination of n = 2 and 3, as written by ; for replacing the a -bT expression. [13] In the present paper, the notation ''HS'' is used for this original method proposed by Hillert and Selleby. Taking the ratio of the second (n = 2) and third (n = 3) terms in this equation as a variable x, the HS method is extended giving the following generalized expression: ensures that the entropy is zero at 0 K. This equation reduces to a -bT at T = T 1 .

Low Temperature Description Below Liquidus
In the third generation of CALPHAD databases, [5,[9][10][11] the temperature dependence of specific heat for pure elements up to the melting point consists of the Einstein function, several power series of temperature and the magnetic contribution in the case where it is applicable: ðEq 7Þ The first term is the contribution from the harmonic lattice vibration, where h E and R represent Einstein temperature and gas constant, respectively. The second term describes the contribution from the electronic excitations and loworder anharmonic corrections. The third term represents contribution from the high-order anharmonic lattice vibrations, where different power series may be used, e.g., i = 2, 3 and/or 4. [2,5,43] The last term is the magnetic contribution, which is zero for pure Al since this element does not show any magnetic ordering. [1,18] Then, the Gibbs energy expression is derived as follows: Applying the expression of Eq 3, the lattice stability is represented as: Note that the values of a 0 and b 0 are different from those of a and b because of the Einstein function. The parameters A, B and C correspond to the Gibbs energy difference attributed to E 0 , to that of the electronic excitations and loworder anharmonic corrections, and to that of the high-order anharmonic lattice vibrations, respectively. All these parameters are optimized using the PARROT module in the Thermo-Calc software [44] to fit the Gibbs energies above room temperatures calculated by a -bT in the existing second generation databases and the E 0 values from ab initio calculations available in existing literature.

High Temperature Description Above Liquidus
In the higher temperature range above the melting point where the solid phase is unstable/metastable, the following two expressions (for C P and G) are used for modelling the solid phase in the third generation of CALPHAD databases. [5,[9][10][11] ðEq 12Þ In the present work we use the same power series as in Eq 12 for the polynomial function in the lattice stability at high temperatures above the melting point, giving the following equation: The parameters I, J, K, L and N are determined to fulfill the conditions that the Gibbs energy, entropy, enthalpy, specific heat and first temperature derivative of specific heat taken from the low-temperature Eq 10 and high-temperature Eq 13 have identical values at the melting temperature. That is the same criteria employed in the development of third generation databases for the unary systems. [5,[9][10][11] 3 Results and Discussion We demonstrate how to apply the Eq 10 to approximate metastable lattice stabilities through a case study of Al, in which the fcc structure is stable up to its melting point (934 K). The SGTE unary description was assessed based on thermodynamic and phase stability of corresponding unary and relevant higher order systems, which indeed forms the basis of the second generation databases. The Gibbs energies of hcp and bcc phases relative to the stable fcc phase at 298.15 K and 10 5 Pa are expressed as [1] : In the following sections, first we show the method for determining the lattice stability of the metastable hcp Al, for which the Debye temperature is known. Secondly, the lattice stability of the (dynamically unstable) bcc phase in Al is shown as a case for which the Debye temperature is undetermined.

Lattice Stability of hcp Al
Contributions from the anharmonic lattice vibration and electronic excitation are rather small and thus the Debye temperature is particularly important at temperatures lower than the Debye temperature h D . The h D can be calculated from experimental or theoretical information of the equilibrium Wigner-Seitz radius r 0 , bulk modulus B and elastic constants (and then Poisson's ratio m), which can be represented as [45][46][47][48][49] : and h, k B and M are reduced Planck constant (Dirac constant), Boltzmann constant and atomic weight, respectively. h D (-3) means the low temperature limit of the Debye temperature. [47] It has been reported that the hcp structure has r 0 = 0.1582 nm, B = 75.02 GPa and the bcc structure has r 0 = 0.1578 nm and B = 78.36 GPa from ab initio calculations. [40] Assuming the same values for m and h D (0)/h D (-3) ratio for the fcc and hcp structures, and using h fcc E = 294 K, [50] the Einstein temperature that is proportional to h D (0), [51] of the hcp Al is estimated to be h hcp E = 288 K.
In the following paragraph, we test how the HS model works for the thermodynamic properties including the effect of the T 1 temperature. To determine the parameters a 0 and b 0 in Eq 4 0 -6 0 , the sum of the first and second terms of Eq 10 are approximated by a linear expression in a specific temperature range as follows: Then, the lattice stability is re-written as Therefore, the corresponding linear expression In the present case of unary Al, the Einstein temperatures, h fcc E = 294 K and h hcp E = 288 K, yield a 00b 00 T = -7.9 -0.547T (298 K \ T \ 2900 K) for the linear approximation in Eq 18. Consequently, we obtain (aa 00 ) -(bb 00 )T = (5481 ? 7.9) -(1.8 -0.547) T = 5488.9 -1.253T, which corresponds to the set of (a 0 = 5488.9, b 0 = 1.253). Figure 1 shows (a) Gibbs energy, (b) entropy, (c) enthalpy and (d) specific heat of the hcp phase relative to the fcc phase in Al, calculated by using the SGTE expression a -bT or using Eq 10 with different T 1 temperatures ranging from 834 to 1234 K and x = 0.5. Strong deviations are observed at high temperatures for all properties and the HS cannot be readily applied to all temperature ranges. It is obvious that the curvatures at high temperatures are sensitive to the selected T 1 value. In addition, it is expected that the ratio of n = 2 and n = 3 terms (i.e., the value of x in Eq 3-6) affects the description of these thermodynamic properties.
The parameters A, B and C in Eq 10, which are functions of T 1 , x, a 0 , and b 0 , were optimized to fit the Gibbs energies calculated by a -bT in the temperature range from 298.15 K to T max (T max : 1934-7034 K). The obtained parameters are listed in Table 1 and the Gibbs energy, entropy, enthalpy and specific heat of the hcp phase relative to the fcc phase are shown in Fig. 2, together with those from a -bT. During the optimization, the Einstein temperature of the hcp Al is fix at h hcp E = 288 K. The HS tells us the relation between the optimized curve and a -bT. In the generalized model Eq 10 and 4 0 -6 0 , the parameters a 0 , T 1 and x can be calculated according to the optimized value of A, B and C (see Table 1), b 0 value being fixed to be 1.253. It is interesting to note that x varies in the range of 0.4 and 0.5 when the Gibbs energy calculated by a -bT is considered between 298.15 K-T max K, which is close to 0.5 suggested in the original HS model, [13] and that a 0 is close to 5488.9. The T 1 temperature increases with increasing T max . The trend of shift in each thermodynamic property against the T 1 temperature is similar between Fig. 1 and 2. As shown in Fig. 2(a), the Gibbs energy decreases with the increase in the T max temperature and the deviation from a -bT becomes larger at low temperatures. It demonstrates how the choice of different T max affects the zero point energy.
The difference in E 0 values between the hcp and fcc structures predicted by ab initio [40,[52][53][54][55] are also shown in Fig. 2(a) and (c). These data differ by up to about 2000 J/mol. If we try to fit the highest value of the ab initio calculation by Boettger and Trickey, [52] it turns out that T max = 2900 K (the green curve) is the optimal choice. However, we prefer to choose E 0 = 3333.7 J/mol obtained by fitting the total energies reported by Mishra [40] to the Murnaghan equation of the state [56] because it is consistent with the prediction of Debye temperature mentioned above.
For the case of hcp Al, if we try to fit the Gibbs energy data calculated by SGTE's a -bT in the full temperature range between 298.15 K and T max K, the resulted E 0 value seems overestimated compared to the ab initio calculations, even the T max temperature is set as high as 7034 K. Conversely, the choice of E 0 = 3333.7 J/mol implies the poor fitting of Gibbs energy to SGTE's a -bT at low temperatures. Therefore, during the optimization the low-temperature Gibbs energies were given a low weight and the parameters were optimized mainly using the Gibbs energies at temperatures between 1527 K (the hcp/bcc transformation temperature in the SGTE database) and 4000 K as well as E 0 = 3333.7 J/mol from the ab initio calculation. [40] The resultant E 0 is close to the predicted value by the ab initio calculation. It is also found that the ratio x is very different from 0.5 when the lattice stability of SGTE is largely overestimated or underestimated. The Gibbs energy curve falls into a range between that with n = 2 and that with n = 3 in Eq 2 when x is between 0 and 1 (see Fig. 1a of Ref 13), but x = 1.214 in the present case means that the curve is slightly lower out of the case with n = 2. Figure 3 shows the Gibbs energy difference between hcp and fcc Al at 0-4000 K. While the deviation from a -bT is not large at high temperatures for T max = 3934, 4934, 7034 and 4000 K, it unreasonably increases at high  temperatures for T max = 1934 and 2900 K. In such cases, another expression, Eq 13, for the high temperature range is necessary. The partitioning of the Gibbs energy description into two temperature ranges is discussed in the following section for the case of bcc Al.
Considering the E 0 value, the expression for T max = 4000 K is recommended in this work. The procedure of determination of lattice stability for hcp Al is summarized as follows: (1) determination of Einstein temperature, (2) optimization of parameters A-C in Eq 10 using the Gibbs energy calculated by SGTE data in a selected temperature range, (3) determination of appropriate temperature range judged by E 0 , and (4) calculation of parameters I, J, K, L and N in Eq 13 for high temperature range, if necessary. a-bT 1934 K 2900 K 3934 K 4934 K 7034 K 4000 K 1996Boe (LDA) [52] 2000Jaf (LDA) [53] 2000Jaf (GGA) [53] 2007Mis (GGA) [40] 2015Lin (GGA) [55] 2004Wan (GGA) [54] Fig. 3 Gibbs energy difference of hcp Al relative to fcc Al at temperatures from 0 to 4000 K calculated by Eq 10 beyond the melting temperature, in which the parameters were optimized using Gibbs energy data from a -bT at 298 \ T \ T max (T max : 1934-7034 K) or at 1527 \ T \ 4000 K denoted as 4000 K

Lattice Stability of bcc Al
This is a case in which the Debye temperature is unknown or unreliable for some reason, such as limited information. The bcc Al is dynamically unstable at 0 K under ambient pressure. [25][26][27][28] Although the Debye temperature is not definable and the entropy has no physical meaning for an unstable phase, one needs to determine the lattice stability in the CALPHAD method. One option is that we define the lattice stability by changing a -bT of the SGTE form using Eq 3 without the Einstein model. Here, we tried expressing the lattice stability of bcc Al using the Einstein model, which is consistent with other structures. An empirical rule for predicting the Debye temperature for a dynamically unstable structure has been proposed, [49] but reasonable result of Gibbs energy at finite temperatures was not obtained for the bcc Al.
If the Poisson ratio of the bcc structure is assumed to be the same as the fcc structure, the Einstein temperature h E is calculated equal to 277 K. The parameters were optimized using the SGTE Gibbs energy a -bT at 298 -2900 K and the result shows that the E 0 is lower than the ab initio calculation results (Fig. 4), implying that the estimated Einstein temperature is too low.
Therefore, the Einstein temperature h E for the bcc Al was optimized as well as other parameters in Eq 10. In the first trial, only the Einstein function was optimized using the SGTE Gibbs energy a -bT at 298-2900 K and E 0 , [40] and the h E became about 243 K. Although the fit is good, as shown in Fig. 4, the difference in anharmonic and electronic excitation is not considered; thus, this h E may be considered as the lowest limit. Thereafter both the Einstein temperature h E and the parameters A, B and C were optimized to fit the ab initio calculated E 0 , [40] the SGTE Gibbs energy by a -bT at 298-2900 K, and the data point of specific heat at 300-900 K from the calculation for h bcc E = 277 K. The optimized h bcc E is 252 K and the E 0 agrees with the ab initio result. [40] Moreover, the difference 1996Boe (LDA) [52] 2000Jaf (LDA) [53] 2000Jaf (GGA) [53] 2007Mis (GGA) [40] 2015Lin (GGA) [55] 2004Wan (GGA) [54] 0  Table 2.
The expression of Eq 10 results in an unreasonable increase in the Gibbs energy at high temperatures, such as 4000 K. In such cases, another expression, such as Eq 13, should be used for the temperature range higher than the melting point of the stable phase. The parameters I, J, K, L and N in Eq 13 for h bcc E = 252 K are listed in Table 3, where these parameters are mathematically determined by the conditions that the Gibbs energy, entropy, enthalpy, specific heat and its first derivative calculated from the low temperature expression Eq 10 and the high temperature expression Eq 13 are equal at the melting point of the stable phase (934 K). The Gibbs energy, entropy, enthalpy and specific heat for h bcc D = 252 K using Eq 10 at lower temperatures and Eq 13 at high temperatures are shown in Fig. 5 in comparison with a -bT. Each curve is continuous at the melting point of the fcc Al (934 K).
The procedure of determination of lattice stability for bcc Al is summarized as follows; (1) optimization of parameters A-C and Einstein temperature in Eq 10 using the Gibbs energy calculated by SGTE data in a selected temperature range, E 0 and C P , and (2) calculation of parameters I, J, K, L and N in Eq 13 for a high temperature range, if necessary.

Lattice Stabilities of All Allotropes of Unary Aluminum
The optimized thermodynamic parameters for the hcp and bcc Al in the present work are summarized in Table 3. Figure 6 shows the Gibbs energy, entropy, enthalpy and specific heat of the hcp and bcc phases relative to the fcc phase in pure Al. The Gibbs energy of the liquid phase [50] is also shown for reference. The lattice stabilities of the hcp and bcc Al recommended in this work are well fitted to those in SGTE [1] at high temperatures, as shown in Fig. 6(a), and deviation at low temperatures was determined considering the ab initio calculations. [40] The entropy demonstrated in Fig. 6(b) proves to be zero at 0 K. The specific heat of the fcc, [50] hcp and bcc phases is shown in Fig. 7. As mentioned in the introduction, the difference in the specific heat between the fcc and hcp or bcc phases is ignored in the SGTE because of the simple a -bT expression of the lattice stability. Conversely, the specific heat of the bcc Al with lower Einstein temperature has higher value and that of the hcp Al is slightly higher compared with that of the fcc Al in the present work. In more physical sense, the contributions from the electronic excitation and the anharmonic lattice vibration (corresponding to T 2 and T 3 terms) should be theoretically considered. A similar way to the SGTE is adopted that the specific heat at high temperatures is close to each other. The difference in the Einstein function is approximated in the temperature range from 298 to 2900 K to be a 00b 00 T = -49.7 The melting point, enthalpy of transformation and entropy of transformation of the present work are compared with those of the SGTE in Table 4. These thermodynamic properties of the present work are close to those of the SGTE database at high temperatures, i.e., for the fcc/ bcc and fcc/hcp transformations, but noticeably revised at low temperatures, particularly for the transformations involving the hcp phase.

Conclusions
A simple method to assess the lattice stability in the framework of the modelling for the third generation CALPHAD databases was proposed. The expression obeys the third law of thermodynamics, meanwhile, taking full advantage of ab initio calculations at 0 K and those existing descriptions in the second generation data for T [ 298 K. The method was demonstrated using the unary Al system as an example, where the lattice stabilities of the metastable hcp and unstable bcc Al were investigated.
In the case where the Einstein temperature information is available or it can be predicted (e.g., hcp Al), only the other parameters, i.e. A-C in Eq 10, in the HS model for extending the linear function are optimized using the Gibbs energies from the existing SGTE description and the results of ab initio calculations. An appropriate temperature range of the Gibbs energy from SGTE should be chosen for optimization so as to fit other available information, such as E 0 predicted by ab initio calculation. In the other case that the Einstein temperature is unknown (e.g., bcc Al), both the Einstein temperature and HS model parameters are optimized using the Gibbs energies from SGTE and E 0 predicted by ab initio calculation.
The ratio x in the HS expression is close to 0.5 when the lattice stability of SGTE is acceptable in a certain range of   Table 3 fcc hcp bcc  Table 3 at temperatures (a) from 0 to 1000 K and (b) from 0 to 4000 K temperature. It is necessary to use a different expression for the Gibbs energy at temperatures higher than the melting point of the stable phase to avoid strong deviation, in which the parameters I, J, K, L and N in Eq 13 can mathematically be determined by the low temperature Eq 10. This generalized method is simple but flexible enough to estimate the lattice stabilities using the existing second generation data down to room temperature, as suggested by SGTE, and ab initio calculations, e.g., total energy at 0 K. The formation energy of stoichiometric compounds and excess energy of solution phases are often expressed by a -bT in the second generation data especially when thermodynamic information is limited. The present method could be readily applied for these energies to avoid any non-zero vibrational entropy at 0 K.