Xiao-Gang Lu
email: xiaogang@thermocalc.se

Qing Chen

A CALPHAD Helmholtz Energy Approach to Calculate Thermodynamic and Thermophysical Properties of Fcc Cu

Keywords: ab initio, thermodynamics user supplied): Debye-Grüneisen model, thermophysical properties, Cu CALPHAD, Helmholtz energy, Debye-Grüneisen model, Thermodynamics, Thermophysical property, Cu

A CALPHAD Helmholtz Energy Approach to Calculate Thermodynamic and Thermophysical Properties of Fcc Cu

The current CALPHAD technique [START_REF] Kaufman | Computer Calculations of Phase Diagrams[END_REF][START_REF] Saunders | CALPHAD (Calculation of Phase Diagrams): A Comprehensive Guide[END_REF][START_REF] Lukas | Computational Thermodynamics: The Calphad Method[END_REF] is based on the Gibbs energy with temperature and pressure as its natural variables. Because most measurements and applications are carried out at atmospheric pressure, the pressure dependence of the Gibbs energy is often omitted in many thermodynamic databases [START_REF] Dinsdale | [END_REF][START_REF]NIST-JANAF Thermochemical Tables, 4th ed[END_REF] and, as a result, volume and the volume-related thermophysical properties cannot be calculated from such databases.

Exceptions were found in a few cases [START_REF]Thermo-Calc Fe alloys/steel database[END_REF][START_REF] Guillermet | [END_REF][START_REF] Saxena | Properties of Earth and Planetary Materials at High Pressure and Temperature[END_REF], where a pressure-dependence term based on the Murnaghan [9], Birch-Murnaghan [10] or other equations of state (EOS) was added to the Gibbs energy expression. Recently, efforts have been reported on assessing thermophysical properties, such as volume, thermal expansivity and bulk modulus [11,12] etc. and on finding an EOS more suitable for high pressure extrapolations [13][14][15].

The CALPHAD Gibbs energy approach requires many fitting parameters since the temperature dependences are treated with polynomials and assessed separately for heat capacity, thermal expansion, and bulk modulus at atmospheric pressure. Furthermore, due to this practice, the obtained parameters may have no physical meanings, and the intrinsic relations among heat capacity, thermal expansion, and bulk modulus are completely lost.

As a consequence, inconsistent results may be obtained for these properties and abnormal behaviors may occur when using them to make high pressure extrapolations, for example a minimum or negative entropy may appear at or over a certain pressure, which is not justifiable theoretically or experimentally [13,14].

While the CALPHAD technique uses Gibbs energy for the sake of convenience to practical constant temperature and pressure applications, theoretical work adopts naturally Helmholtz energy because ab initio calculations are always performed at constant temperature and volume [see for example, [16][START_REF] Anderson | Equations of state of solids for geophysics and ceramic science[END_REF][START_REF] Herper | [END_REF][19][20]. When the temperature and volume dependence of Helmholtz energy is established with some physical models, other thermodynamic properties including thermophysical properties can be derived in an intrinsically consistent way, and thus abnormal behaviors can be avoided for properties at high temperatures and high pressures. By using this approach, a systematic work has 3 been published for volume and thermal expansivity for transition cubic metallic elements [21], cubic carbides and nitrides [22].

For clarity, we shall now distinguish between two categories of physical properties in this paper. One is the thermodynamic properties, i.e. the properties taking into account in the traditional CALPHAD approach, such as heat capacity, enthalpy, entropy and Gibbs energy. The other is the thermophysical properties, including volume, thermal expansivity, bulk modulus, shear modulus, Young's modulus, and Poisson's ratio.

In this paper, we study element metal Cu and propose to adopt the theoretical Helmholtz energy approach within the framework of CALPHAD technique so that experimental data on both thermodynamic and thermophysical properties can be assessed at the same time, and a set of parameters having physical meanings can be obtained to describe most of experimental information. If no experimental data is available, especially in the case for most metastable phases, ab initio calculation results at 0K can be used directly in this approach to estimate the thermodynamic and thermophysical properties over a certain range of temperatures and pressures. Firstly we introduce a Helmholtz energy approach based on the Debye-Grüneisen model and the free electron Fermi gas model, and then explain the assessment procedure by using the optimization module PARROT in Thermo-Calc [23][24][START_REF] Jansson | TRITA-MAC 0234[END_REF]. In Section 4, we present in detail the assessment results of thermodynamic and thermophysical properties of fcc Cu. A discussion about the facilitated link between CALPHAD and ab initio calculations is given in Section 5, followed by a summary in the end.

Methodology

Different from the traditional CALPHAD Gibbs energy approach [START_REF] Kaufman | Computer Calculations of Phase Diagrams[END_REF][START_REF] Saunders | CALPHAD (Calculation of Phase Diagrams): A Comprehensive Guide[END_REF][START_REF] Lukas | Computational Thermodynamics: The Calphad Method[END_REF], the present approach is based on Helmholtz energy. For a system at constant temperature and volume, its total Helmholtz energy F comprises the total energy at 0 K (i.e. static lattice energy),

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) V T TS V T E V T TS V T E V E V T F V T F V E V T F el el D D tot el D tot , , , , ) , ( , , - + -
+ = + + = (1)
where T is temperature and V is volume. E D and S D are the vibrational energy and entropy, respectively. E el and S el are the counterparts due to electronic excitations. The contributions summed up in the above equation are adequate to accurately describe simple metals like Cu, which will be justified in Section 4. For transition metals, the anharmonicity effect due to phonon-phonon interactions, as well as the magnetic and vacancy contributions may become important and proper models should be devised, which will be considered in subsequent papers. Now we give expressions for the contributions considered in Eq. ( 1) one by one.

Total static energy at 0 K, E tot

The total energy at 0 K is a function of volume or a linear separation (such as lattice parameter which is frequently used in ab initio calculations). It is also called as equation of state (EOS) at 0 K. In this work, it is represented by a Morse function [see, e.g. 16]:

( ) ( ) ( ) ref x x x x c E e e E x E + - - = - - - - 0 0 2 2 ) ( ϕ ϕ , (2)
where x represents a linear separation, usually lattice parameter, x 0 is the equilibrium x, φ is a fitting parameter. In ab initio calculations, E ref is usually chosen as the energy of a system with same number of isolated atoms (where x is infinitely large), and E c is then the cohesive energy since E c =E ref -E 0 at x=x 0 . In this work, in order to be consistent with SER (Stable Element Reference), the standard reference state adopted in CALPHAD community, E ref is chosen to ensure that the enthalpy of the stable phase of an element is zero at 298.15 K and 1 bar. E c is thus no longer equal to the value of the cohesive energy. Subsequent work is needed for the comparison.

2.2 Helmholtz energy of electron gas, F el (T,V)

The thermal excitations of electrons are not negligible for many materials, especially at cryogenic temperatures and at high temperatures near melting points. The simplest way to model the electronic excitations is the free electron Fermi gas model, which treats the conduction electrons as a freely moving electron gas and applies Pauli exclusion principle and Fermi-Dirac distribution.

However, the interaction between lattice vibrations and conduction electrons (i.e. the electron-phonon coupling) in a metal causes a deviation of the electronic thermal behavior from the free electron model. Grimvall [27] derived a temperature dependent function to describe this effect. According to White [28], the electronic contribution at high temperatures can not be extrapolated from the low-temperature linearly even when the electron-phonon interaction disappears at high temperatures, because of the complexity of the electron energy band structure, particularly for the transition metals.

For the case of Cu, we will see in Section 4.3 that the free electron model is adequate to reproduce the available data for fcc Cu. The work by Greeff et al. [29] gives a further proof which will also be shown below. Therefore in this work we will not pursue to model the electron-phonon coupling which will be considered for the transition metals in future work.

For a system of N electrons, the heat capacity of the electron gas is linearly related to T [START_REF] Kittel | Introduction to solid state physics[END_REF],

3 (

) ( 2 2 2 2 T k D T T k N C B F F B el V ε π π = = , (3)
T C el el V γ = , (4)
) 5
In order to relate F el to volume, γ el is treated as a function of volume, but temperatureindependent, i.e.

( )

a el el V V V         = 0 0 γ γ , ( 6
)
where 0 el γ is the Sommerfeld parameter at V 0 , the equilibrium volume at 0 K and 1 bar.

The parameter 'a' is a constant. Both parameters can be evaluated from the measured heat capacity at cryogenic temperatures and 1 bar, and the parameter 'a' is further adjustable to fit high temperature data.

In addition, theoretic evaluations are feasible when D(ε F ) is known. By using Eq. ( 3) and

(4), we can relate el γ to D(ε F ) that is calculated by the following two methods. According to the free electron gas theory, D(ε F ) is expressed as follows, ( )

2 / 1 2 / 3 2 2 2 2 F e F m V D ε π ε       = h , (7) 3
/ 2 2 2 3 2         = V N m e F π ε h , (8)
where m e is the electron rest mass, ћ is the Planck constant. Another way to obtain D(ε F ) is through ab initio electronic calculations, which give more realistic results. Moruzzi et al. [16] evaluated the Grüneisen parameter from ab initio calculated energy and volume data at 0 K to describe the volume dependence of the Debye temperature.

Wang et al. [42] later on improved and extended this method, and made calculations at high temperatures and pressures. Our calculation method was inspired by Wang et al.'s method, and they were shown to be equivalent [21].

There are three well-known approximations to account for the Grüneisen parameter γ, and can be combined into the following expression [42][START_REF] Zharkov | Equations of State for Solids at High Pressures and Temperatures[END_REF][START_REF] Wasserman | [END_REF],

[ ] [ ] V PV V PV V V ∂ ∂ ∂ ∂ - - = + + / / 2 ) 1 ( 3
1 ) ( ) 1 )( 3 / 2 ( 2 ) 1 )( 3 / 2 ( 2 λ λ λ γ , ( 9
)
where

V V E V P ∂ ∂ - = ) ( ) (
. When λ = -1,0 and +1, one obtains the Slater approximation [START_REF] Slater | Introduction to Chemical Physics[END_REF], ( )

2 / 1 3 / 2 ) ( 3 ) 1 ( 2 ) (       + - ∂ ∂ - = V V P V V P DV V D λ θ , ( 10
)
where D is the integral constant that is determined below. Note that P(V) and E(V) concerns only EOS at 0 K.

The Debye temperature can also be related to the sound velocity, υ D [see, e.g. [START_REF] Grimvall | Thermophysical Properties of Materials[END_REF][START_REF] Chen | [END_REF]:

( ) D a B D V N k V υ π θ 3 / 1 2 6         = h , (11)
where N a is the number of atoms and

m VB v k B v k D ) ( ) ( = = ρ υ , ( 12
)
where ν is Poisson's ratio, B the bulk modulus, ρ the density, m the mass, and k(ν) a coefficient depending on ν:

3 / 1 2 / 3 2 / 3 ) 1 ( 3 1 3 1 ) 2 1 ( 3 ) 1 ( 2 3 2 ) ( -                 - + +       - + = ν ν ν ν v k (13)
Inserting the definition of B, one can get:

( ) ( ) 2 / 1 3 / 2 3 / 1 2 2 / 1 2 3 / 1 2 ) ( 6 ) ( ) ( 6 ) (       ∂ ∂ - =           ∂ ∂ -         = V V P V N m k v k m V V P V V N k v k V a B a B D π π θ h h (14)
By comparing the above equation with the Slater approximation of Eq. ( 10), i.e. when λ = -1:

F o r P e e r R e v i e w O n l y 9 ( ) 2 / 1 3 / 2 ) (       ∂ ∂ - = V V P DV V D θ , ( 15
)
we can identify the constant D as:

( ) 3 / 1 2 6 ) ( a B N m k v k D π h = . ( 16
)
Inserting it back into Eq. ( 10), we obtain

( ) ( ) 2 / 1 3 / 2 3 / 1 2 ) ( 3 ) 1 ( 2 ) ( 6 ) (       + - ∂ ∂ - = V V P V V P V N m k v k V a B D λ π θ h . ( 17
)
When V increases to a critical value, the quantity in the bracket may become negative, and θ D is not defined. However, the thermal expansion of normal materials is not big enough to reach this critical value. Furthermore, when the vibrational contribution is taken into account, another case of undefined θ D will occur as discussed in Section 4.9. It also deserves mentioning that in the vicinity of the equilibrium volume, the magnitude of the corresponding pressure is much smaller than that of bulk modulus (hundreds GPa).

Therefore, one can safely neglect the second term in the bracket, i.e. θ D becomes λindependent, when not too high pressure is concerned.

Introducing the general relation, Eq. ( 17), into the energy and entropy expressions of the Debye approximation as follows [START_REF] Grimvall | Thermophysical Properties of Materials[END_REF], one can thus express the vibrational energy and entropy as a function of volume and temperature, ( )

      + = T TD k N k N V T E D B a D B a D θ θ 3 8 9 , , (18) ( )
( )       - -       = - T D B a D D e T D k N V T S / 1 ln 3 4 3 , θ θ , ( 19
)
where D(θ D /T) denotes the Debye function.

The model discussed above takes into account the volume dependence of the Debye temperature, which is within the quasiharmonic approximation and the anharmonicity due to thermal expansion is accounted for as the calculated thermal expansivity is non- ). Wallace [50] showed that for 16 nearly-free-electron elements, the anharmonicity arising from the phonon-phonon interactions is negligible. Especially for Cu, it is nearly zero at the melting temperature. For Cr, Mo and W the anharmonic contributions to entropy and energy are higher and explicit modeling of anharmonicity is needed [50][51][52]. However, this is beyond the scope of the present work.

Calculation of thermodynamic and thermophysical properties

It is known that the pressure is the negative value of the derivative of the Helmholtz energy with respect to volume. Given certain external conditions, i.e. temperature T and pressure P, the equilibrium volume and energy are determined numerically from the Helmholtz energy curve (Fig. 1) by fulfilling the external conditions. The Gibbs energy is then known by G=F+PV. Other thermodynamic and thermophysical properties can in turn be determined by various derivatives of the Helmholtz energy.

Among various properties, the bulk modulus B is calculated from the derivative of the pressure with respect to volume. Poisson's ratio is evaluated from experimental data as shown below. Once the bulk modulus and Poisson's ratio are known, for polycrystal or homogeneous isotropic materials we can calculate Young's modulus, E, and shear modulus, S, by:

) 1 ( 2 ) 2 1 ( 3 ν ν + - = B S , ( 20
) ) 2 1 ( 3 ) 1 ( 2 ν ν - = + = B S E . ( 21
)
3. Assessment procedure PARROT is a powerful tool in Thermo-Calc to optimize model parameters [23]. In principle, all types of experimentally determined thermodynamic and thermophysical data at various temperatures and pressures can be considered into PARROT and optimized at the same time. In practice, some experimental data are excluded in the should also be adjusted manually during the optimization. The parameters x 0 , E c and φ in Eq.2 are the key parameters optimized in PARROT using all selected experimental data.

The parameter E ref in Eq.2 is a constant and trivial if the SER reference is not observed for a single phase. For a multiphase system, it has to be assessed to account for the relative phase stability. Finally, all experimental data are used to fine-tune the model parameters. The minimum corresponds to the equilibrium free energy and volume at P=0. As the temperature increases, the minimum shifts towards higher volume and lower Helmholtz energy. The static lattice refers to the lattice at 0 K without any zero-point lattice vibrations. When the zero-point vibrations are considered, small positive contributions to both the total energy and equilibrium volume occur. The latter gains about 0.6% at 0 K for fcc Cu.

Heat capacity at 1 bar, C p

In Fig. 2, the calculated heat capacity at 1 bar is compared with the experimental data from Brooks et al. [53], Martin [54], Stevens and Boerio-Goates [55], as well as the data compiled by Dinsdale [START_REF] Dinsdale | [END_REF] and the NIST-JANAF data [START_REF]NIST-JANAF Thermochemical Tables, 4th ed[END_REF]. The NIST-JANAF data adopted the selected data by Hultgren et al. [START_REF] Hultgren | Selected values of the Thermodynamic properties of the elements[END_REF] In the metastable region above the melting temperature, the NIST-JANAF data are generally higher than the present prediction. The electronic contribution at high temperatures is also shown in Fig. 2. The effect is small but noticeable.

The electronic contribution to heat capacity

At sufficiently low temperatures, the Debye T 3 law is valid for the vibrational heat capacity C D , i.e.

Adding the electronic contribution from Eq. ( 4), we obtain a linear relation,

2 / / AT T C T C el V P + = ≅ γ , ( 23
)
where A is a material related constant, and γ el is the Sommerfeld parameter and treated as volume-dependent and temperature-independent (see Eqs. ( 4) and ( 6)). To obtain γ el , two parameters, i.e. 0 el γ and 'a' in Eq. ( 6), are assessed from the measured heat capacity data at cryogenic temperatures [55,[START_REF] Osborne | [END_REF][58][59] as well as theoretical evaluations. In the course of assessments, we firstly obtain electronic density of state (DOS) at Fermi energy by ab initio calculations. For a set of volumes in the vicinity of the equilibrium volume, the electronic DOS's were calculated by Lu et al. [21]. The DOS data are used in Eqs. ( 3) and ( 4) to evaluate γ el as a function of volume. A recent ab intio calculation by Greeff et al. [29] presents a similar result (if a suspect misprint of the unit is corrected).

In the second step, the obtained γ el and experimentally measured C p and V m data are used to assess 0 el γ and 'a' in Eq. ( 6). The value of 0 el γ obtained, i.e. 0.695 mJ/mol/K 2 , is the same as the value obtained by Phillips [62] and later listed by Kittel [START_REF] Kittel | Introduction to solid state physics[END_REF]. The result can be read from the intersection with the vertical axis in Fig. 3. The results for γ el are shown in Fig. 4.

In Section 2.2, we discussed the electron-phonon coupling but commented that the free electron model is able to reproduce the experimental data for fcc Cu. Miiller and Brockhouse [63] calculated the total lattice heat capacity (harmonic plus anharmonic) at constant pressure from the measured frequency -wave vector dispersion relations for the lattice vibrations. The electronic C p up to 900 K was then derived by taking the difference between the total lattice C p and the experimentally determined C p . They concluded that the derived electronic C p agreed well with the linear relation (Eq. 4) for temperatures lower than 700K, above which the derived data were higher than the linear relation. It is worth noting that Eq.4 describes the linear behavior of C v and the corresponding C p should be higher. In Fig. 5, similar results are obtained but we adopt recent assessed data by Dinsdale [START_REF] Dinsdale | [END_REF] and NIST-JANAF data [START_REF]NIST-JANAF Thermochemical Tables, 4th ed[END_REF] for the experimentally observed C p in the calculations. The agreement is acceptable. Also in the figure, our previously reported data [21] calculated by the following equations are shown as dotted line.

(

( ) ∫ ∫ - = F d V n N d f V n N V T E el ε ε ε ε ε ε ε ε 0 ) , ( ) , ( , )
( ) ( ) ( ) [ ] ∫ - - + - = ε ε d f f f f V n Nk V T S B el 1 ln 1 ln ) , ( , , (24)
Where n(ε,V) is the electronic DOS, and f(ε) represents the Fermi-Dirac distribution. [79], as well as the AIP [START_REF] Gray | American Institute of Physics Handbook[END_REF] and TPRC data [START_REF] Touloukian | Thermal Expansion of Metallic Elements and Alloys[END_REF]. Almost the same accuracy is obtained as the previous assessment by Lu et al. [12], and the values can now be reproduced down to 0 K.

Molar volume and the coefficient of linear thermal expansion (CTE

Elastic properties for polycrystal at 1 bar

The adiabatic bulk modulus for polycrystal is experimentally available through measuring the adiabatic single crystal elastic constants c 11 , c 12 and c 44 , which are then estimated by the VRH (Voigt-Reuss-Hill) approximation [START_REF] Grimvall | Thermophysical Properties of Materials[END_REF]. Alternatively the bulk modulus is derived from the measurement of sound velocity. Due to the presence of crystal defects, the experimental data for elastic moduli are scattered as shown in Fig. 8 for bulk modulus. Ledbetter [80] measured the sound velocities for polycrystalline copper at 295K, and compared with the average value from 18 single-crystal experiments. The adiabatic bulk modulus calculated from his sound velocities is 138.9 GPa, while the VRH approximation gives 138.1 GPa. These values were assigned higher weights than other data at room temperature during the assessment. Our calculated value, 137.9 GPa, agrees well with Ledbetter's data, while the data from other sources [81][82][83][84] below room temperature are slightly lower. At high temperatures, the calculated adiabatic bulk modulus is in good agreement with the measured data by Chang and Himmel [85].

According to Ledbetter, the calculated adiabatic shear moduli at 295 K are 47.86 and 46.87GPa, respectively. In the present approach, the shear modulus is calculated from the assessed adiabatic bulk modulus and Poisson's ratio according to Eq. (20 In the work by Gerlich and Kennedy [86], the shear modulus from zero pressure to 2.0GPa was reported. However, no temperature is indicated in the paper and the measured data at zero pressure were about 10% higher than Ledbetter's data and our result.

Properties at high pressures

The molar volume of fcc Cu at high pressures and room temperature was measured by Dewaele et al. [87] and the pressure was calibrated by two standards. The present calculated result lies between these two calibrations as plotted in Fig. 9. In addition, the normalized volume agrees well between the present calculation and the measurements [88][89][90] as shown in Fig. 10.

The pressure derivative of the isothermal bulk modulus B' at 298K calculated in the present work is 5.10, which is lower than that derived from the measurement by Klooster 12, and the shear modulus is calculated along the Hugoniot using Eq. ( 20) as shown in Fig. 13.

As already mentioned above, the conventional CALPHAD approach to model highpressure properties always causes abnormal property behaviors at high temperatures or high pressures. Using the present calculation approach, all calculated properties behave reasonably. Two examples are presented up to 1000GPa in Fig. 14. Note that at 1000 GPa the quasiharmonic approximation and the present approach may already fail but Fig. 14 is only for the purpose of demonstration. In the present approach, we do not try to represent the complexity but calculate the Debye temperature at finite temperatures from bulk modulus (or sound velocity), Poisson's ratio and volume by:

( ) D a B D V N k T υ π θ 3 / 1 2 6         = h , (26)
where

m VB v k D ) ( = υ
. Note that Eq. ( 26) is formally the same as Eq. ( 11), a formula which emphasizes the volume (instead of temperature) dependence of θ D and involves only the properties at 0 K. In Eq. ( 26), however, the temperature dependence of θ D is considered through the temperature dependence of bulk modulus and volume. We emphasize that the Poisson's ratio in both Eq. ( 11) and ( 26) to calculate D υ is a constant value, i.e. 0.358 for fcc Cu. This type of temperature variation of θ D is merely an approximation of the complex temperature-dependence of Debye temperature, but adequate to describe most of the properties as we have shown above. One exception is the heat capacity between ~ 20 K and 200 K, where the agreement is obviously not as good as at high temperatures (see Fig. 2). In the following discussion we show the possibility to examine heat capacity, the corresponding 'heat capacity θ D ' at low temperatures. By examining Eq. ( 26) and knowing that the bulk modulus and volume are determined by experimental data, it is convenient to treat Poisson's ratio (i.e. the parameter ν given in Table 1) as a fitting parameter when calculating θ D . By adjusting ν and fitting the lowtemperature Cp, the 'heat capacity θ D ' is calculated and plotted in Fig. 16, together with the results evaluated by Holste et al. [58]. The imposing character of the 'heat capacity θ D ' is a sharp decrease at cryogenic temperature before recovery. On the other hand, a value of 0.358 is assigned to Poisson's ratio (see Table 1), giving a best overall agreement with all types of experimental data. The corresponding θ D is plotted as a dashed curve in Fig. 16. Between ~ 20K and 200K, the dashed curve is above the 'heat capacity θ D ', which is exactly the cause of the imperfection of the agreement of the heat capacity at low temperatures. In fact in the inset in Fig. 2, the calculated heat capacity is slightly lower than the measured data, implying that the calculated θ D is higher. One can

Grüneisen parameter

The thermal Grüneisen parameter, γ G , is experimentally observable through measuring the adiabatic decompression (∆T/∆P) S :

P S S S G C VB P T T B α γ 3 =       ∆ ∆ = , ( 27
)
where α is the coefficient of linear thermal expansion (CTE), other symbols have their usual meanings.

Ramakrishnan et al.

[96] measured (∆T/∆P) S for Cu at high pressures, and calculated γ G from Eq. ( 27) before fitting it as a function of volume and extrapolating to very large compression. In Fig. 17, we compare our calculations with the measured (∆T/∆P) S . In Fig. 18, γ G is calculated in a large pressure range and plotted as a function of compression by using Eq. ( 27).

High-temperature instability

In Fig. 1, the Helmholtz energy is plotted against volume, where minima (corresponding to P=0) exist on the curves for temperature below 2230 K. Above 2230K, the Helmholtz energy decreases monotonously with increasing volume and no minimum can be found.

This means the bulk modulus is negative and the Debye temperature is not defined. The direct consequence is that the heat capacity, thermal expansivity as well as other properties can not be calculated above 2230 K. In this way, we can avoid the situations that the heat capacity becomes infinitely large at the high-temperature metastable region and solids reappear above their melting points. This type of 'instability' is different from the usually discussed 'mechanic instability', where the shear elastic modulus C' defined al. [21,22] have applied this method to study some metallic elements including Cu and transition-metal carbides and nitrides.

Cares should be taken due to the fact that many of the 'metastable' states necessary for the CALPHAD modeling are actually unstable, with undefined entropies resulting from strong instabilities due to phonon softening, etc. To justify the predictions, one should always take into account ab initio calculated elastic constants and full phonon spectra, or check Bain deformation paths.

For stable phases with plenty of measurements, the present calculation approach allows a direct comparison between the assessed and ab initio calculated EOS at 0 K. In Fig. 19, the energy is normalized by subtracting the minimum energy E 0 , and the lattice parameter for the fcc lattice is subtracted by the lattice parameter a 0 at E 0 . The reason is that in ab initio calculations by Lu et al. [21], the reference energy is different from the present assessment, and ab initio calculated lattice parameter usually deviates slightly from the observation. As seen from Fig. 19 for fcc Cu, the agreement is excellent in view of the fact that the assessment does not use the ab initio data. With the growing interests in extending thermodynamic databases to high-pressure region, as well as in modeling thermophysical properties, the CALPHAD methodology faces challenges. Fundamental modifications are proposed to model thermodynamic and thermophysical properties in an internally consistent manner within the framework of the Helmholtz energy. This approach has long been adopted in theoretical researches.

The proposed calculation method pulls closer the CALPHAD assessments and ab initio calculations. Since the main parameters in the models are physically sound, their magnitudes can be estimated from ab initio calculations and then assessed in the PARROT module. On the other hand, when experimental information is not available, using ab initio calculated EOS at 0 K as an input, this approach predicts valuable information at finite temperatures and pressures (while cares should be taken by further studying instability issues).

Although the Helmholtz energy approach is successfully applied to fcc Cu in the present work, the difficulty in extending it to multi-component multi-phase alloying systems should not be underestimated. It is not clear at this stage how difficult it is to apply this method to complex alloy systems, and how one can account for the magnetic contribution and the ordering phenomena. It is also desirable to develop a suitable model for liquids.

Having said that, the proposed approach has obviously a number of advantages:

1) The thermodynamic and thermophysical properties, such as the heat capacity, entropy, volume, thermal expansivity and elastic moduli, are modeled consistently by using only a few parameters. The traditional CALPHAD modeling does not couple those properties together.

2) No abnormal behavior is observed for the properties at high temperatures and pressures.

3) The model parameters possess physical meanings. The assessments can thus benefit greatly from ab initio calculations and can be compared directly with theoretical work.

4) It is possible to account for the low-temperature regime down to 0 K. Fig. 1 Helmholtz energy varies with volume at different temperatures for fcc Cu. The dotted curve is the energy for a static lattice excluding zero-point energy, while the solid curve at 0K corresponds to the vibrating lattice including the zero-point contribution. The inset is a zoom in the curves at around the minimum energy at 0 K and 300 K. Fig. 2 The calculated heat capacity at atmospheric pressure compared with the experimental data for fcc Cu. The imperfection shown in the inset is explained in the text. Fig. 3 The calculated heat capacity of fcc Cu at cryogenic temperatures. The solid and dashed lines are the calculation results with and without electronic contributions, respectively. The intersection of the solid line with the vertical axis gives 0 el γ in Eq.6.

F

Fig. 4 The Sommerfeld parameter for fcc Cu is determined as a function of volume. The dashed line and the triangles are evaluated from the electronic DOS obtained from ab initio calculations, and the dotted line is evaluated from the electronic DOS in the free electron model. Fig. 5 The calculated electronic heat capacity at atmospheric pressure compared with the derived experimental data for fcc Cu (see text for details). Squares: Dinsdale's data [START_REF] Dinsdale | [END_REF] were used for total Cp. Triangles: NIST-JANAF data [START_REF]NIST-JANAF Thermochemical Tables, 4th ed[END_REF] were used. Fig. 19 The EOS for energy (E) and lattice parameter (a) at T=0K. E and a are normalized by the minimum energy E 0 and corresponding a 0 , respectively. The two end points on the calculated curve (solid curve) correspond to P= -20 GPa (right) and 200 GPa (left), while at the minimum energy P=0.

4 E

 4 tot , the contribution of the vibrating lattice, F D , and the contribution due to the thermal excitations of electrons, F el :

  to compare the Morse function with other frequently used equation of state functions, such as Murnaghan[9], Birch-Murnaghan[10] and Vinet function[START_REF] Vinet | [END_REF].

  and using Eq. (9), one obtains the dependence of the Debye temperature D θ on volume V:

4. 1

 1 The total Helmholtz energy In Fig.1, the Helmholtz energy as a function of volume at various temperatures is plotted.

  from 300K to the melting point. The assessment is in good agreement with the NIST-JANAF data below the melting temperature, and with the recent data from 200 K to 400 K by Martin, and by Stevens and Boerio-Goates.Close to the melting temperature, the compiled data by Dinsdale are slightly lower than the present calculation and other data sources. If the data by Brooks et al. were fitted, the calculated thermal expansivity would be much higher than the experimental data. Under the present approach, Brooks et al.'s data are not favored. In fact, the NIST-JANAF data did notice Brooks' data were 1-2% higher. At low temperatures shown in the inset of Fig.2, the agreement is not perfect. The reason is explained in Section 4.7 where the Debye temperature and Poisson's ratio are discussed.

  et al. [84], i.e. 5.33. By solely fitting Dewaele et al.'s data using the conventional ruby calibration, the calculated B' is 5.25. A bigger B' requires that the calculated volumes at high pressures fit with Dewaele et al.'s data by the new calibration method other than the ruby calibration. Probably the present calculated volumes at very high pressures and B' are slightly low. The Hugoniot shock wave curve, with 298K and 1 bar as initial conditions, is calculated. The temperature and volume along the Hugoniot are shown in Fig.11. In the figure for volume, the onset of melting on the Hugoniot is marked at about 230 GPa according to Mitchell et al. [91-92]. Since in the present work liquid Cu is not considered, the experimental Hugoniot data above 230 GPa can not be reproduced. The shock temperature was calculated by Peng et al. [93] using the Grüneisen parameter, specific heat and density, and agrees well with the present calculation. From the measured sound velocity in shock-loaded copper, Hayes et al. [94] estimated the shear modulus and Poisson's ratio among other elastic properties. The zero pressure ratio is 0.356, which is in a good agreement with the present assessed value. Peng et al. [93] calculated the shear modulus below 100 GPa using the same sound velocity data as Hayes et al.. In the present work, Poisson's ratio is fitted to the data by Hayes et al. as shown in Fig.

4. 7

 7 Debye temperature, θ D , and Poisson's ratio, ν The Debye model approximates the complex phonon frequency spectra with a parabolic relation between the phonon density of state and lattice frequency. Fig. 15 demonstrates the comparison between the Debye phonon spectra with θ D =320 K, and the measured phonon frequency spectra by Larose and Brockhouse [95]. The Debye cutoff frequency ω D is determined by θ D . Clearly only for long wave length (i.e. low frequency) the Debye spectra satisfactorily represent the experimental phonon spectra, while ω D is close to the experimental maximum phonon frequency. In addition, the phonon spectra vary with temperature.Therefore, the common applications of the Debye model in the literature, using a singlevalued θ D to represent experimental data over a wide temperature range, are usually not satisfactory. To appropriately reproduce the experimental data at different temperatures, θ D should change with temperature in a complex way. Furthermore, one often refers to, for example, 'entropy θ D ', 'heat capacity θ D ' et al. when studying data on entropy and heat capacity, respectively[START_REF] Grimvall | Thermophysical Properties of Materials[END_REF].

  behavior of Cp below 20 K, but it is not prominent because of the dominant electronic contribution at the cryogenic temperatures.

  Figure captions

Fig 6

 6 Fig 6 Calculated molar volume of fcc Cu at atmospheric pressure compared with the experimental data.

Fig. 7

 7 Fig.7The calculated coefficient of linear thermal expansion (CTE) at atmospheric pressure compared with the experimental data for fcc Cu.

Fig. 8

 8 Fig.8The calculated adiabatic (Bs) and isothermal (B T ) bulk modulus at atmospheric pressure compared with the experimental data for fcc Cu.

Fig. 9

 9 Fig.9The calculated molar volume of fcc Cu at 298 K and high pressures.

Fig. 10

 10 Fig. 10 Normalized molar volume of fcc Cu at 298K at high pressures.

Fig. 11

 11 Fig.11The temperature (a) and volume (b) along the Hugoniot. The dashed line marks the onset of melting at about 230 GPa[92], above which concerns about liquid that is not considered in the present work.

Fig. 12

 12 Fig.12Poisson's ratio along the Hugoniot. Note that it is not calculated from the present approach but fitted to Hayes et al. data.

Fig. 13

 13 Fig.13 Shear modulus along the Hugoniot.

30 Fig. 14

 3014 Fig.14 The calculated isobaric heat capacity and coefficient of linear thermal expansion (CTE) at high temeratures and high pressures. No abnormal behavior is found.

Fig. 15

 15 Fig.15 The measured phonon frequency spectra for fcc Cu at 296 K, compared with the Debye phonon spectra with θ D = 320 K.

Fig. 16

 16 Fig. 16 Debye temperature calculated by adjusting ν and fitting experimental heat capacity data (see Section 4.7 for details).

Fig. 17

 17 Fig. 17Calculated and measured adiabatic (∆T/∆P) S . This value is measurable and used to evaluate the Grüneisen parameter.

  Fig. 17Calculated and measured adiabatic (∆T/∆P) S . This value is measurable and used to evaluate the Grüneisen parameter.

Fig. 18 Calculated

 18 Fig.[START_REF] Herper | [END_REF] Calculated Grüneisen parameter at different volumes (pressures) at room temperature (solid line). The triangles are evaluated from the measured (∆T/∆P) S in Fig.17using Eq.27, while the dashed line is the extrapolation from the measured data by Ramakrishnan et al. [96].

  is the Fermi temperature, k B is Boltzmann constant, D(ε F ) is the electronic density of state (DOS) at the Fermi energy ε F . Equivalently, we can use the so-called Sommerfeld parameter, γ el , to express the heat capacity and Helmholtz energy,

	where T F
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	6

  on compared with the calculations to further confirm the optimization results. In the present optimization, we use the experimental data on the heat capacity, volume (at 1 bar and high pressures), thermal expansivity and bulk modulus. In case many measurements are available for the same type of property, personal judgments based on experimental accuracy are required in order to select suitable experimental data to fit all data reasonably or identify inconsistent data.All the parameters needed to be assessed are listed in Table1, together with the values for fcc Cu. The two parameters for the electronic contributions, i.e. 0

	optimization but later
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	11

el γ and 'a' defined in Eq.6, are firstly evaluated from experiments at cryogenic temperatures and ab initio calculations. Since at cryogenic temperatures the electronic contribution is dominant, a set of reasonable initial values for other parameters is adequate for a good assessment of the two parameters. Secondly, the model for the Grüneisen parameter, i.e. Slater's model or DM model, must be chosen before assessing other parameters. The corresponding parameter, λ, is thus not a fitting parameter in the optimization procedure. Only three options are allowed for the Grüneisen parameter model, i.e. λ = -1, 0 and +1. One has to manually fix it by trial and error. It is worth mentioning that the last option, which corresponds to the free volume theory, may not be a good choice in most cases since it leads to very low values for thermophysical properties. Preferably, the Poisson's ratio ν,

Table . 1

 . Model parameters that can be assessed or input and their values for fcc Cu.

	φ		1.73097	-	Fitting parameter in Morse function (Eq.2),
					optimized by PARROT.
	E ref	2.892033E+05 J/mol	Reference energy to maintain SER (Stable
					Element	Reference)	equivalency,	i.e.
					H 298K,1bar =0 ( see Eq.2). Assessed separately
					after fixing x 0 , E c and φ.
	ν		0.358	-	Poisson's ratio, fixed and adjusted manually.
	λ		0	-	Approximation	options	for	esimating
					Grüneisen parameter: Slate: -1; DM: 0 (See
	γ a	0 el	F 6.95E-4 o r 0.60	Eq.9) J/mol/K 2 Sommerfeld parameter at 0K and 1 bar (see Eq.6), assessed separately using Cp at cryogenic temperatures. -For describing volume dependence of
			P e r e	Sommerfeld parameter (see Eq.6), determined from ab intio calculations, and fine-tuned by using experimental C p and V m data.
				R e v i e w
						O n
							l y
	Parameter for fcc Cu	Unit	Note
	x 0	3.5956E-010	m	Lattice parameter in Morse function (Eq.2),
					optimized by PARROT.
	E c	2.972229E+05 J/mol	Equivalent cohesive energy in Morse function
					(Eq.2), optimized by PARROT.
					12

4. Results and discussions for fcc Cu

Cu is probably one of the most thoroughly studied elements. Both experiments and ab initio calculations are accumulated to allow us to carry out a detailed study by using the proposed approach in this work.

All calculated properties for fcc Cu are shown in Figs.

1-18

, along with the corresponding experimental data. A summary of the thermodynamic and thermophysical properties at 298.15 K and 101325 Pa is given in Table

2

. The parameters for fcc Cu assessed in this work are listed in Table

1

.

Table . 2

 . Thermodynamic and thermophysical properties at 298.15K and 101325 Pa calculated in the present work. lower than experimental data[80]. If the Poisson's ratio is 0.345, the shear modulus and Young's modulus are 47.68 and 128.25 GPa, respectively.

	Coefficient of linear thermal expansion,	16.52 10 -6
	α, K -1
	Adiabatic bulk modulus, B s , GPa	137.90
	Isothermal bulk modulus, B T , GPa Adiabatic shear modulus, S s , GPa Adiabatic Young's modulus, E s , GPa	133.97 43.26 (a) 117.49 (a)
	Poisson's ratio, ν	0.358
	Heat capacity Debye temperature, θ D , K 319.13
	Sound velocity, m/s	2440.72
	o Grüneisen parameter, γ G F (a)	1.99
	r
	P
	e e r
	R e v i e w
		O n
		l y
	Isobaric heat capacity, C p , J/mol/K	24.46
	Isochoric heat capacity, C v , J/mol/K	23.76
	Enthalpy, H m , J/mol	0.0
	Entropy, S m , J/mol/K Molar volume, V m , m 3 /mol Lattice parameter, a, m	32.34 7.1103 10 -6 3.6146 10 -10
	Density, ρ, kg/m 3	8937.2
		13

  Osborne et al.[START_REF] Osborne | [END_REF] presented the experimental heat capacity from 1 to 25 K, and compared with 15 measurements in other laboratories. Holste et al.[58], Hurley and Gerstein[59] later presented a result in agreement with Osborne et al.'s data. Tsumura et al.[60] studied the heat capacity below 1 K. Swenson[61] measured the heat capacity of Cu below 30 K and found that the agreement with Holste et al. data is extremely good.

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	However, no experimental value was presented in his paper. The calculated Cp in this
	work is plotted in Fig.3.
	15

  The calculated molar volume and CTE at 1 bar are shown in Figs.6 and7, respectively, and compared with the experimental data [64-75 for volume, 76-79 for CTE]. The best agreement for CTE is obtained with the data provided by White and Minges

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	) at 1 bar
	16

  =1/2(c 11 -c 12 ) is negative. The latter is irrelevant to fcc Cu, but for metastable or unstable states it is an important issue to study.

	as C'5. Application of ab initio calculations to the CALPHAD method
	Ab initio calculations become more and more involved in the CALPHAD modeling, from
	studying fundamental problems like lattice stability to providing physical properties to F assist database development; from 0 K ground-state calculations to finite-temperature o vibrational excitation research. Using the present approach in PARROT, a new research
	r method is provided to easily fit the ab initio calculated energy-volume (or lattice
	parameter) EOS at 0K, and effectively predict the thermodynamic and thermophysical P properties at finite temperatures and high pressures. These data are especially valuable e r e for metastable phases, or when experimental data are limited at extreme conditions. Lu et
	R e v i e w
	O n
	l y
	21

Table . 1

 . Model parameters that can be assessed or input and their values for fcc Cu.

	[93] J.Peng, F.Jing, D.Li, J Appl Phys 98(2005) p.013508.
	[94] D.Hayes, R.S.Hixson, R.G.McQueen, Shock Compression of Condensed Matter-
			1999, edited by M.D.Furnish, L.C.Chhabildas, R.S.Hixson (American Institute of
			Physics, New York, 2000), p.483.
	[95] A.Larose, B.N.Brockhouse, Can J Phys 54(1976) p.1990.
	[96] J.Ramakrishnan, R.Boehler, G.H.Higgins, G.C.Kennedy, J Geophys Res 83(1978)
			o F o F o o F p.3535.
			r r r r
			P P P P
	e e r e e r e e r Parameter for fcc Cu Unit e x 0 3.5956E-010 m e r E c 2.972229E+05 J/mol	Note Lattice parameter in Morse function (Eq.2), optimized by PARROT. Equivalent cohesive energy in Morse function
	φ E ref ν λ	R e v i e w R e v i e w R e v i e w (Eq.2), optimized by PARROT. Fitting parameter in Morse function (Eq.2), R optimized by PARROT. 2.892033E+05 J/mol 1.73097 -Reference energy to maintain SER (Stable e Element Reference) equivalency, i.e. v H 298K,1bar =0 ( see Eq.2). Assessed separately i after fixing x 0 , E c and φ. e 0.358 -Poisson's ratio, fixed and adjusted manually. 0 -Approximation options for esimating w Grüneisen parameter: Slate: -1; DM: 0 (See
	γ a	0 el	6.95E-4 0.60	O n l y O n l y O n l y O J/mol/K 2 Sommerfeld parameter at 0K and 1 bar (see Eq.9) Eq.6), assessed separately using Cp at cryogenic temperatures. n -For describing volume dependence of l Sommerfeld parameter (see Eq.6), determined y from ab intio calculations, and fine-tuned by
					using experimental C p and V m data.
					27

Table . 2

 . Thermodynamic and thermophysical properties at 298.15K and 101325 Pa calculated in the present work.

	Isobaric heat capacity, C p , J/mol/K	24.46
	Isochoric heat capacity, C v , J/mol/K	23.76
	Enthalpy, H m , J/mol	0.0
	Entropy, S m , J/mol/K Molar volume, V m , m 3 /mol Lattice parameter, a, m	32.34 7.1103 10 -6 3.6146 10 -10
	Density, ρ, kg/m 3	8937.2
	Coefficient of linear thermal expansion,	16.52 10 -6
	α, K -1
	Adiabatic bulk modulus, B s , GPa F o Isothermal bulk modulus, B T , GPa Adiabatic shear modulus, S s , GPa Adiabatic Young's modulus, E s , GPa	137.90 133.97 43.26 (a) 117.49 (a)
	r Poisson's ratio, ν	0.358
	Heat capacity Debye temperature, θ D , K 319.13 Sound velocity, m/s 2440.72 P Grüneisen parameter, γ G 1.99 (a) lower than experimental data [80]. If the Poisson's ratio e is 0.345, the shear modulus and Young's modulus are 47.68 e r and 128.25 GPa, respectively.
	R e v i e w
		O n
		l y
		28

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters