Anharmonic Properties of Aluminum from Direct Free Energy Interpolation Method∗

Zhi-Guo Zhao(赵志国),Jun-Sheng Sun(孙俊生),Xiu-Lu Zhang(张修路),Hai-Feng Yang(杨海峰), and Zhong-Li Liu(刘中利),*

1College of Physics and Electric Information,Luoyang Normal University,Luoyang 471022,China

2The Unit 63615 of People’s Liberation Army,Kuerle 841001,China

3Laboratory for Extreme Conditions Matter Properties,Southwest University of Science and Technology,Mianyang 621010,China

1 Introduction

Like platinum and tungsten,aluminum is another frequently used pressure standard material in high-pressure experiments,such as the shock wave experiment.The pressure-volume-temperature equation of states(EOS)provides key information for it as a pressure standard.[1]Al is a typicalsp-bonded simple metal and crystallized in face-centered-cubic structure under ambient conditions.With the simple electronic and lattice structures,it is often taken as a typical prototype for theoretical calculations.[2−3]For the theoretical calculation of free energy,the accurate anharmonic effects are very hard to be taken into account in the theoretical EOS in previous methods.

Usually,the quasi-harmonic approximation(QHA)was used to calculate the temperature effects and it only takes account into the volume dependence of phonon frequencies of lattice vibrations.However,in real materials the phonon frequencies depend largely on temperature,especially in the highly anharmonic solids or at the temperature close to melting.This is because that the QHA omits part of anharmonicity and is only accurate when temperature below Debye temperature.[4−6]In order to overcome these diffculties of QHA,several anharmonic free energy calculation techniques have been developed,such as the thermodynamic integration technique[7−8]and the self-consistent ab initio lattice dynamics(SCAILD)method.[9−10]

Molecular dynamics(MD)simulations include anharmonicity to any order beyond Debye temperature,[11−12]especially applicable to high-temperature free energy calculations.Thus it can correct the free energy of the QHA considerably in the free energy calculations.In this paper,we apply the MD simulations to extract anharmonic free energy of Al and then construct more accurate EOS beyond the QHA.

The rest of the paper is organized as follows.In Sec.2,we describe the details about the computation method.We illustrate the results in detail and make deep discussions in Sec.3.The conclusions are presented in Sec.4.

2 Computational Method

2.1 Quasi-Harmonic Approximation

The free energy of Al was first calculated by the quasi-harmonic approximation and then corrected by taking into account the phonon-phonon interactions.The quasi-harmonic approximation(QHA)describes volumedependent thermal effects based on lattice vibration models known as phonon,neglecting the phonon-phonon interactions.In the framework of the QHA,the Helmholtz free energy of a crystal system is,

whereEstatic(V)is the zero-temperature energy of a staticlattice with volumeV.The termFzpis the zero-point motion energy written as,

whereg(ω,0K)is the 0 K phonon density of states(PDOS)of the frequencyω.

The termFph(V,T)in Eq.(1)is the lattice vibrational free energy,and is calculated from,

2.2 Anharmonicity Corrections

In the harmonic approximation,the PDOSg(ω,0K)is calculated from the force constants using either the direct method[13−14]or the density functional perturbation theory.[5,15]In the harmonic approximation,the phononphonon interaction contribution to free energy are neglected,resulting in large errors of thermodynamic properties at high temperature.Therefore,the phonons interactions must be fully taken into account to include the rest of the anharmonic effects.However,the rest anharmonic free energy calculation is very challenging.As we know that the explicit anharmonicity are naturally implied in molecular dynamics(MD)simulations.We can take into account the rest anharmonic effects from MD simulations.The PDOS atT,g(ω,T),can be extracted from the Fourier transform of autocorrelation of atomic position function r(t)in a constant temperature MD simulation.[16]

According to the Wiener–Khintchine theorem,[17−18]the autocorrelation of position function r(t)is given by

The PDOS atTis,

The Helmholtz free energy of a crystal system is then rewritten as,

whereEstatic(V)is the total energy of a static lattice with volumeV.

In Eq.(6),the termFzpis the zero-point motion energy of the lattice written as,

The termFph(V)in Eq.(6)is calculated from,

After obtaining the finite number of temperature dependent PDOS(TD-PDOS)data,we use a much straitforward and simple method named the direct free energy interpolation(DFEI)method,to construct full anharmonic effects at any temperature.

The DFEI method extracts anharmonic free energy at any temperature based on the limited several TD-PDOS data from MD simulations.The free energiesF(V,T)can be directly deduced from the TD-PDOS data at the specific temperaturesTs with the temperature interval 50 K,according to Eq.(6).The free energies at any temperatures were obtained by interpolating these resulting free energy data using the numerical analysis technique of the basis spline(B-spline).Finally,the anharmonic effects of other thermodynamic properties were corrected.

2.3 The Interatomic Potential of Al

We perform classical MD simulations to simulate Al at different pressure and temperature.Then,we gather the trajectories of all the atoms in crystal to calculate phonon density of states.The embedded-atom-method(EAM)[19−20]potential developed by Mishin et al.[21−22]was used to describe the interatomic interactions of Al atoms.For a metal system containingNatoms,the total potential energyEtotis the sum of the embedding energyFand a pair potentialϕ,

whererijis the distance between the atomsiandj.The functionFi(ρi)is the energy to embed the atomiinto the background electron densityρi,which is the superposition of the atomic densities,

TheMD simulationswereconducted with the large-scale atomic/molecular massively parallel simulator(LAMMPS)[23−24]package.The simulation box constructed from the multiplication 14×14×14 of the facecentered-cubic(fcc)conventional unit cell including 10976 atoms.The simulations were conducted for all the supercells with different volumes in the canonical ensemble(NVT).The periodic boundary condition was used for all the atoms in the simulation box.The time step was 1 fs and the total number of time steps were 10000.From Fig.1,we see that 10000 MD steps are able to converge the PDOS well.

Fig.1 (Color online)The PDOS obtained using different numbers of MD steps.

3 Results and Discussions

3.1 High-Temperature Phonon Dispersion and DOS

The phonon dispersion curve of Al were calculated using the direct method from molecular dynamics simulations.[25]The calculated phonon dispersion curve of Al are compared well with experimental data at 298 K[26]in Fig.2(a).The frequency of lattice vibration decreases with temperature until very small negative frequencies occur at 900 K as shown in the circle of Fig.2(a).Such very small negative frequencies often result from the numerical noises in the direct phonon calculation method with molecular dynamics simulations.[26]The good agreement of the calculated phonon dispersion curve with experimental data show the validity of the EAM potential for Al.The TD-PDOSs of Al with atomic volumeV0=16.61˚A3are obtained by Fourier transform of autocorrelation of atomic position function and plotted in Fig.2(b).The PDOSs are from 0 to 900 K with the interval of 100 K.The density of states of higher frequencies increase with temperature,attributed to the increasing collective excitations of higher frequencies from the phonon-phonon interactions.On the contrary,the lower frequencies of vibrations are inhibited with increasing temperature.The 900 K PDOS does not show any negative frequency and thus does not affect the free energy calculation.

The volumes of supercells used for PDOS calculations had volumes varying from 1.08V0to 0.64V0with the interval of 0.1V0,whereV0is the equilibrium volume at ambient pressure.The simulation temperature for each volume ranged from 0 K to 950 K with the interval of 50 K.The volume range and interval made a good sample for the description of the equation of states of Al(Fig.3(a))compared with the experimental data,[27]indicating that the anharmonic free energy interpolation based on the PDOS with the temperature interval of 50 K is suffcient.

Fig.2 (a)The comparison of the calculated 0 GPa phonon dispersion curves of Al with experimental data at 298 K.[26]The circle indicates the negative frequencies at Γ point.(b)The phonon DOS of Al at atomic volume of 16.61˚A3from 0 to 900 K with the interval of 100 K along the arrow direction.

3.2 Thermal Expansion Coeffcients

We calculate the Helmholtz free energy from the PDOS according to Eq.(1).The third-order Birch–Murnaghan equation of states[28]was fitted to the calculated free energyF–Vdata at each temperature. The third-order Birch–Murnaghan equation of states is given by,

whereη=(V0/V)1/3,andF0,V0,B0andB′are fitting parameters.The thermodynamic properties and the anharmonic effects are automatically analyzed using the our PhaseGO package.[29−30]

The accuracy of the calculated thermal properties was checked by the thermal expansion coeffcient defined by,The calculated thermal expansion coeffcients as a function of temperature are shown in Fig.3(b).The QHA results are in good agreement with experimental data[31]below 400 K,while when the temperature is above 400 K the zero-pressure thermal expansion coeffcient gradually deviates from experiment.This attributes to the neglect of anharmonicity caused by phonon-phonon interactions.This deviation is also observed in the calculations for other transition metals,such as Ta,[32]Pt,[33]and Pd.[34]However,the DFEI results show good agreement with the experimental data[31]in Fig.3(b).

Fig.3 (a)The comparison of the calculated isothermal compression curve with experimental data.[27](b)The calculated thermal expansion coeffcients of Al compared with experimental data.[31]

3.3 Heat Capacity

Fig.4 The specific heat of Al at constant volume and at constant pressure compared with the experimental data.[35]

whereFis the Helmholtz free energy calculated from Eq.(6).The thermal expansion caused by anharmonic effects results in a difference betweenCPandCV.The difference between the two can be written as

whereαVis the volume thermal expansion coeffcient andB0is the bulk modulus.Because the heat capacities are

The specific heat at constant volume was calculated by obtained from a second derivative of the Helmholtz free energy,they are more sensitive to the tiny errors in the Helmholtz free energy.

Figure 4 showsCVandCPas a function of temperature at 0 GPa.The calculatedCVaccords reasonably well with experimental data[35]below 600 K,and first increases dramatically as pressure increases and then finally approaches to 3R.The DFEI and QHA results are very similar.The calculatedCPis also in good comparison with experimental data.[35]The better agreement with experimental data is also found after anharmonic correction(DFEI),especially whenTis beyond 800 K.

3.4 Grüneisen Parameter and Bulk Modulus

The thermodynamic Grüneisen parameter was achieved by

whereBTis isothermal bulk modulus.The isothermal bulk modulus was calculated from

The calculated isothermal bulk modulus versus temperature curves at different pressures are shown in Fig.5(a).The isothermal bulk modulus curve at 0 GPa is very good agreement with the experimental data.[36]The calculated thermodynamic Grüneisen parameter is compared with the ab initio data from Ref.[37]in Fig.5(b).The calculated thermodynamic Grüneisen parameter is also in very good agreement with the reference data from Ref.[38].

Fig.5(a)The calculated isothermal bulk modulus of Al at different pressures,in comparison with zero-pressure data[36](open circles).(b)The calculated thermodynamic Grüneisen parameter compared with ab initio data from Ref.[37].

4 Conclusion

In conclusion,we calculate and correct the thermodynamic properties of Al using the direct calculation method of full free energy of lattice vibrations,i.e.,the DFEI method which makes accurate anharmonic corrections for the thermodynamic properties of materials beyond the QHA.The calculated thermal expansion coeffcients of Al using DFEI method agree better with experiment than the QHA results.Furthermore,the calculated isotherms,the constant volume and constant pressure heat capacities,the bulk moduli using DFEI are all in good agreement with experimental data.This indicates that the DFEI method is indeed an effcient and accurate method to determining the thermodynamic properties of materials from the high temperature phonon density of states.

[1]C.J.Pickard and R.J.Needs,Nat.Mater.9(2010)624.

[2]R.Martin,Nature(London)400(1999)117.

[3]J.C.Boettger and S.B.Trickey,Phys.Rev.B 53(1996)3007.

[4]U.Hansen,P.Vogl,and V.Fiorentini,Phys.Rev.B 60(1999)5055.

[5]S.Baroni,S.D.Gironcoli,A.D.Corso,and P.Giannozzi,Rev.Mod.Phys.73(2001)515.

[6]S.Baroni,P.Giannozzi,and E.Isaev,Rev.Miner.Geochem.71(2010)39.

[7]D.Alf`e,G.D.Price,and M.J.Gillan,Phys.Rev.B 64(2001)45123.

[8]D.Alf`e,G.D.Price,and M.J.Gillan,Phys.Rev.B 65(2002)165118.

[9]P.Souvatzis,O.Eriksson,M.I.Katsnelson,and S.P.Rudin,Phys.Rev.Lett.100(2008)95901.

[10]P.Souvatzis,O.Eriksson,M.I.Katsnelson,and S.P.Rudin,Comput.Mater.Sci.44(2009)888.

[11]I.Errea,M.Calandra,and F.Mauri,Phys.Rev.B 89(2014)64302.

[12]R.Car and M.Parrinello,Phys.Rev.Lett.55(1985)2471.

[13]D.Alf`e,Comput.Phys.Commun.180(2009)2622.

[14]A.Togo,F.Oba,and I.Tanaka,Phys.Rev.B 78(2008)134106.

[15]X.Gonze,Phys.Rev.B 55(1997)10337.

[16]M.Thomas,M.Brehm,R.Fligg,et al.,Phys.Chem.Chem.Phys.15(2013)6608.

[17]N.Wiener,Acta Math.55(1930)117.

[18]A.Khintchine,Math.Ann.109(1934)604.

[19]M.S.Daw and M.I.Baskes,Phys.Rev.B 29(1984)6443.

[20]S.M.Foiles,M.I.Baskes,and M.S.Daw,Phys.Rev.B 33(1986)7983.

[21]J.M.Winey,A.Kubota,and Y.M.Gupta,Modell.Simul.Mater.Sci.Eng.17(2009)55004.

[22]J.M.Winey,A.Kubota,and Y.M.Gupta,Modell.Simul.Mater.Sci.Eng.18(2010)29801.

[23]S.J.Plimpton,J.Comput.Phys.117(1995)1.

[24]http://lammps.sandia.gov.

[25]L.T.Kong,Comput.Phys.Commun.182(2011)2201.

[26]R.Stedman and G.Nilsson,Phys.Rev.145(1966)492.

[27]A.Hnstrm and P.Lazor,J.Alloy Compd.305(2000)209.

[28]F.Birch,Phys.Rev.71(1947)809.

[29]Z.L.Liu,Comput.Phys.Commun.191(2015)150.

[30]Z.L.Liu,Comput.Phys.Commun.197(2015)341.

[31]F.C.Nix and D.MacNair,Phys.Rev.60(1941)597.

[32]Z.L.Liu,L.C.Cai,X.R.Chen,et al.,J.Phys.:Condens.Matter 21(2009)95408.

[33]T.Sun,K.Umemoto,Z.Wu,J.C.Zheng,and R.M.Wentzcovitch,Phys.Rev.B 78(2008)24304.

[34]Z.L.Liu,J.H.Yang,L.C.Cai,F.Q.Jing,and D.Alf`e,Phys.Rev.B 83(2011)144113.

[35]M.Forsblom,N.Sandberg,and G.Grimvall,Phys.Rev.B 69(2004)165106.

[36]A.Dewaele,P.Loubeyre,and M.Mezouar,Phys.Rev.B 70(2004)94112.

[37]J.L.Nie,L.Ao,F.A.Zhao,M.Jiang,and X.T.Zu,Can.J.Phys.93(2015)55004.

[38]W.B.Holzapfel,M.Hartwig,and W.Sievers,J.Phys.Chem.Ref.Data 30(2001)515.