Liquid Aluminum: atomic diffusion and viscosity from ab initio molecular dynamics.

We present a study of dynamic properties of liquid aluminum using density-functional theory within the local-density (LDA) and generalized gradient (GGA) approximations. We determine the temperature dependence of the self-diffusion coefficient as well the viscosity using direct methods. Comparisons with experimental data favor the LDA approximation to compute dynamic properties of liquid aluminum. We show that the GGA approximation induce more important backscattering effects due to an enhancement of the icosahedral short range order (ISRO) that impact directly dynamic properties like the self-diffusion coefficient. All these results are then used to test the Stokes-Einstein relation and the universal scaling law relating the diffusion coefficient and the excess entropy of a liquid.

Increasing demands to optimize methods and facilities for material processing technologies such as melting, refining, or casting metallic alloys are boosting the precise knowledge of thermophysical properties of these alloys. The main goals are an improvement of the final product quality, an enhancement of the process efficiency, and an economical consumption of resources and energy. More particularly, as the solidification process of a liquid alloy has a strong impact on the structure and properties of the solid material, the understanding of the behavior and more particularly the knowledge of the thermophysical properties of molten alloys prior to solidification are essential for the development of materials with predetermined characteristics. For instance, it is known that properties of lightweight aluminum-based alloys are affected by the temperature conditions of melting and casting1, and the knowledge of the diffusion process in their liquid phases is hence a necessity.

Among all thermophysical properties, diffusion coefficients are known to play a prominent role for a detailed understanding of the solidification process including crystal growth234 and vitrification5. Quite surprisingly, experimental values for liquid aluminum and its alloys are scarce. Their direct measurements in liquid alloys are mostly made using so-called capillary tube methods that require isotopes as tracers. However, the influences of convective flow on the diffusion profile during annealing can be a big problem6. As a consequence, self-diffusion coefficients D can be overestimated and published values often differ by a factor of 2. Moreover, radioactive isotopes are not available for aluminum. D are also often obtained indirectly from experimental determinations of the shear viscosity coefficients, η, easier to measure and the use of the Stokes-Einstein (SE) relation7, DSE = kBT/2π Rη, kB being Boltzmann's constant and R is an effective atomic diameter taken as the position of the first peak of the pair-correlation function. If the applicability of the SE relation to liquid metals is well accepted above their melting temperatures, the scatter of experimental viscosity values is important in the case of aluminum8, reflecting difficulty of obtaining accurate data. This can be due to the fact that Al has a high affinity to oxygen. Solid oxide particles can possibly be trapped between the melt and the crucible, affecting drastically the flow profile and therefore the viscosity. For instance, the scatter between measurements made with different techniques is ±14% in the case of aluminum while it is only ±6% for iron8.

Recently, it has been shown for some liquid metals like Ni, Ti, or Cu91011 that self-diffusion coefficients can be obtained accurately from the use of quasielastic neutron scattering (QNS) in a large temperature range, via the determination of the incoherent intermediate scattering function. For liquid Al, the situation is more complicated since its incoherent cross section is two orders of magnitude lower than that of liquid Cu. Very recently, Demmel et al.12 and Kargl et al.13 report temperature-dependent self-diffusion coefficients of liquid Al, with different assumptions to extract the incoherent contribution from the scattering function. The two sets of experimental data are in good agreement and display an Arrhenius-type behavior. The absolute values of D at and around the melting temperature agree also well with those obtained from the Stokes-Einstein relation. However, at high temperatures, i.e. 1200 K, the discrepancy becomes important12.

From a theoretical point of view, molecular dynamics (MD) is a widely used simulation technique for investigating the structural and dynamic properties in the liquid state14. In classical molecular dynamics, the simulated particles are moved using forces extracted from semi-empirical interatomic potentials, electronic degrees of freedom remaining hidden. By contrast, ab initio molecular dynamics (AIMD) within the density functional theory (DFT) are based on quantum-many-body forces computed via the Hellmann-Feynman theorem, which is straightforward using plane-wave basis15.

Properties of liquid Al have been studied extensively using the semi-empirical embedded-atom method (EAM) formalism. However Becker et al.16 have shown that diffusion coefficients are strongly dependent on the different implementations of the EAM potential. Moreover, the best agreement leads to an underestimation of the experimental values that becomes more important at high temperatures.

Within the DFT framework, Alfe and Gillan17 have shown that the calculated self-diffusion and shear viscosity coefficients are sensitive to the choice of the density. Their calculations gave D values within the range 0.52–0.68 Å2/ps while η values are in the range 1.4–2.2 mPa.s. Gonzalez et al.18 used an orbital-free DFT method to compute dynamic properties of liquid Al at T = 943 and 1323 K. This approach underestimates the experimental data, the discrepancy becoming very important at high temperatures. Such a result shows that the accuracy of the orbital-free DFT approach to compute dynamic properties of liquid Al is likely to be limited by the use of an approximate electronic kinetic energy functional and a local pseudopotential. In a more recent work, the same group19 performed Kohn-Sham ab initio molecular dynamics simulations and reported the same diffusion coefficient D as in Ref. 17. They also computed the shear viscosity from the transverse current correlation function. They obtained η = 1.05 mPa.s, which is smaller than the value obtained by Alfe and Gillan17, using the definition of η as the time integral of the auto-correlation function of the off-diagonal components of the stress tensor.

In this paper, we address the important issue of the determination of the dynamic properties of liquid Al using ab initio molecular dynamics simulations. Both the local-density approximation2021 (LDA) and the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof22 (PBE) form for the exchange-correlation energy were investigated. We show that the determination of the self-diffusion coefficient as well as the viscosity depends on the exchange-correlation functional used in DFT calculations. More particularly, GGA based calculations of the diffusion coefficient underestimates both LDA based values and experimental data. We demonstrate that this behavior is related to the detailed description of the short-range order and more particularly to the enhancement of the icosahedral short-range order using the GGA approximation. All these results are then used to test the Stokes-Einstein relation and the universal scaling law relating the diffusion coefficient and the excess entropy of a liquid2324.

Results

In this section, we determine the single atom as well collective dynamics that give access to transport properties such as the self-diffusion coefficient, D, and viscosity, η.

Self-diffusion coefficient

In a first step, we consider the single-atom dynamics through the velocity auto-correlation function, which is defined as: where ν(t) are the position and velocity of atom i at time t, and the angular brackets denote the average over time origins as well as atoms. The time integration of the velocity auto-correlation function gives access to the self-diffusion coefficient, D. This is equivalent to the determination of D using the linear slope at long times of the mean-square displacement725, and require shorter simulation times to have statistically meaningful results.

The temperature dependence of D using the LDA and GGA approximations is reported in Fig. 1. We can see that GGA calculations underestimate LDA calculations by about 20% above the melting point and by about 30% in the undercooled state. LDA calculations are closer to the experimental values1213 and similar to previous LDA-AIMD simulations at T = 1000 K1719. Assuming an Arrhenius-type behavior for the diffusion process in the liquid state, an activation energy E can be derived. As shown in Fig. 1, the fit describes the data quite well in the whole range of temperature. The derived activation energy is E = 249 ± 5 and 260 ± 5 meV respectively for the LDA and GGA approximations. Activation energies have been also obtained from experimental data. The reported values are 280 ± 70 (Ref. 13) and 274 ± 4 meV (Ref. 12) and are only slightly larger than our calculated values.

Figure 1

Evolution of the self-diffusion coefficients of liquid Al as a function of inverse temperature.

The solid circles are AIMD results with LDA and the open ones are those with GGA. The dashed lines are their respective Arrhenius fit. The open triangle and the star correspond to previous ab initio calculations with LDA911. The open diamonds and full squares are experimental data of Refs. 6 and 7, respectively, along with their Arrhenius fit (dash-dotted lines).

Viscosity

Another interesting dynamical magnitude is the transverse current correlation function C(q,t) which has the advantage of yielding a generalized q-dependent shear viscosity from which the hydrodynamic limit can be evaluated2627. It is a straightforward approach based on the decay of the correlation of the atom motions, which is directly linked to the viscosity, and does not require the use of the forces as needed in the Green-Kubo formalism28. We briefly describe here the main points of the formalism. Firstly, the function C(q,t) is defined as where J(q,t) represents the transverse current that is expressed along the x direction as The quantity ν,(t) is the x-component of the velocity of atom i and q is a wave vector along the z direction, recalling that q = (2π/L)(n, n, n) are reciprocal vectors which have to be compatible with the length L = V1/3 of the simulation cell, and n, n and n are integers. Two formally identical expressions are also used for y and z directions. The zero limit Laplace transform gives rise to the q-dependent shear viscosity k and m stand respectively for Boltzmann's constant and mass of the atoms. In the hydrodynamic limit, i.e. by extrapolating η(q) to q = 0, the shear viscosity η can be obtained. The extrapolation is performed with the function proposed by Alley and Alder26 which is used in general to represent the viscosity of a dense hard-sphere packing and has proven to work well for liquid Ni29. Fig. 2 shows that the fitting gives a good representation of the LDA and GGA results of q-dependent viscosity for liquid Al at different temperatures.

Figure 2

q-dependent viscosity of liquid Al obtained from our AIMD with (a) LDA, and (b) GGA.

The symbols are the calculations using the transverse current correlation function for T = 875 K or 850 K (open diamonds), T = 1000 K (open squares), T = 1125 K (open triangles), and T = 1250 K (open circles). The solid lines are the fits using Eq. (5).

At 1000 K, the LDA calculated value is 1.17 mPa.s close to the assessed value2 of 1.18 mPa.s. The values given by previous AIMD simulations using LDA are 1.4 and 1.05 mPa.s, respectively obtained from the time integral of the auto-correlation function of the off-diagonal components of the stress tensor17 and from the transverse current correlation function19. The value afforded by the Stokes-Einstein relation is 1.21 mPa.s, by using R = 2.72 Å, which corresponds to the first peak of the LDA pair-correlation function. The good agreement between this value and the value of 1.17 mPa.s calculated above supports the validity of using the SE relation to describe dynamic properties of liquid Al using LDA. For GGA, the value obtained from the SE relation is 1.53 mPa.s (R = 2.72 Å corresponding to the first peak of the GGA pair-correlation function) which can be compared to that obtained by the direct method, namely 1.37 mPa.s. Both values overestimate the assessed value and their difference is larger than that obtained in LDA. This also holds for the other temperatures, as displayed in Fig. 3. It can be seen that the results of LDA are close to the assessed values of Ref. 8. Moreover, we note that they are also close to the experimental values of Iida and Shirashi30 and remain consistent with those of Battezzati and Greer31, both experimental data sets being retained in the compilation by Dinsdale and Quested32.

Figure 3

Viscosity of liquid Al as a function of temperature.

The full circles and full squares correspond to our AIMD results, respectively with LDA and GGA, as obtained with the transverse current correlation functions. The open circles and open squares correspond to viscosities inferred from the self-diffusion coefficients, respectively with LDA and GGA, determined from our AIMD results using the Stokes-Einstein relation. The open triangle and the star correspond to previous ab initio calculations911. The dashed, dash-dotted, and dash-dot-dotted lines correspond to the experimental values of Ref. 23 and 24, and the assessment of Ref. 2, respectively.

Discussion

Differences between GGA and LDA calculations of the self-diffusion coefficient and the viscosity are significant, more than 15% for the viscosity and 20% for diffusion. To understand these discrepancies, we inspect some computed structural properties using both GGA and LDA functionals.

In Fig. 4, we display the pair distribution function g(r) of liquid Al for all the thermodynamic states as described in the preceding section. We can see that our simulations agree well with experimental X-ray diffraction measurements3334 as well as with previous ab initio simulations19 for the liquid state. We do not see any appreciable difference between the two sets of calculations but it is known that changes which may occur in the short ranger order (SRO) require more detailed analysis3536.

Figure 4

Evolution of the pair-correlation function with temperature.

The open circles and triangles correspond to experimental values of Refs. 26 (1023 K) and 27 (1123 K and 1223 K), respectively. The curves for 1125 K, 1000 K and 875 K are shifted upwards by an amount of 1, 2 and 3, respectively.

Integration of the computed pair distribution function g(r) up to its first minimum gives the coordination number, which is displayed in Fig. 5(a) as function of temperature. We find an increase of coordination number from N = 11.3 ± 0.1 at T = 1250 K to N = 11.8 ± 0.1 at T = 875 K while GGA calculations give N = 11.6 ± 0.1 at T = 1250 K and N = 12.1 ± 0.1 at 850 K. We obtain that LDA calculations give coordination numbers smaller than those obtained using GGA calculations but differences remain small and cannot explain differences obtained in the dynamic properties. Note that the experimental values of N reported by Mauro et al.34 increase from 12.9 ± 0.2 to 13.1 ± 0.2 from T = 1123 K to T = 1273 K. However the authors of Ref. 33 estimated that these values are overestimated by about 8%. Correction of the experimental coordination numbers by this amount would lead to values in good agreement with our AIMD results and with a similar temperature evolution as shown in Fig. 5(a).

Figure 5

(a) Evolution of the coordination number with temperature. The solid squares are the results of LDA and the open ones are those of GGA. The open triangles are the corrected experimental values of Ref. 18; (b) Evolution of abundances of the most important pairs for liquid Al with temperature. The solid symbols are the results of LDA and the open ones are those of GGA. Error bars are typically of the order of 0.01.

More insight into the structural changes is gained by analyzing the inherent structures and characterizing the local environment surrounding each atomic pair that contributes to the first peak of g(r) in terms of the number and properties of common nearest neighbors of the pair under consideration35. In Fig. 5(b), we report the most abundant pairs, averaged over the ten inherent configurations for each temperature and for both LDA and GGA calculation, i.e., 142× (sum of 1422 and 1421), 1431 and 15×× (sum of 1551 and 1541) pairs found in both calculations. The number of 15×× bonded pairs is a direct measure of the degree of icosahedral ordering including both perfect and distorted icosahedral motifs while 142× bonded pairs are characteristic of close packed structures. The 1431 pairs either can be considered as distorted icosahedra or distorted close-packed structures37383940. As they are similar in both approximations and hold steady with temperature, they are not considered to be responsible for differences between dynamic properties using either GGA or LDA. The same comment holds for the 142× pairs. On the contrary, Fig. 5(b) shows that the number of 15×× pairs is more important using the GGA approximation. Our findings indicate that the icosahedral short range order is more pronounced when using the GGA approximation. We can notice also that the icosahedral 15×× fraction always exceeds the close-packed 142× fraction and it is the only fraction that displays a linear increase as the temperature decreases. This trend is in fine agreement with results based on the atomistic cluster alignment method41.

To correlate ISRO with dynamic properties, we focus on the backscattering regime that is predominant for pure liquid Al because of its high density. The backflow induced by a moving atom increases the probability of an atom to jump back toward its initial position. It is characterized by a well-pronounced minimum in the velocity auto-correlation function which is very sensitive to the local liquid density. In Fig. 6, we report the velocity auto-correlation function computed at T = 1000 K using the LDA and GGA approximations. The first minimum becomes deeper using GGA, which reveals an increase of backscattering effects. Note that this difference becomes more important when the temperature decreases, as shown in the inset. This variation can be related to the increase of the icosahedral motifs using the GGA approximation, as discussed in Ref. 42. Icosahedral motifs are known to be the most compact local structures and consequently the backscattering is more pronounced in the liquid which presents the highest ISRO and the diffusivity is smaller.

Figure 6

Velocity auto-correlation function obtained for liquid Al at 1000 K with LDA (solid lines) and GGA (dashed lines).

Inset: same comparison in the vicinity of the main minimum for 875 or 850 K, 1000 K, 1125 K, and 1250 K. For the last three temperatures the curves are shifted upwards by an amount of 0.05, 0.1 and 0.15, respectively.

Understanding the connection between the transport coefficients and structural properties of dense liquids can be also tackled using the correspondence between the transport coefficients and the internal entropy established by Rosenfeld24 or Dzugutov23. The latter author proposed a universal scaling relationship between the excess entropy of a liquid, S, and a dimensionless form of the diffusion coefficient, D*. The excess entropy is the total entropy of the liquid minus the ideal gas contribution. The relationship can be written as: Where D* is given by:In Eq. (7), D is the diffusion coefficient and Γ is the Enskog collision frequency which can be calculated for the temperature T and the density ρ as: where m and σ are the atomic mass and the hard-sphere diameter. g(σ) represents the value of the pair correlation function g(r) at the contact distance.

In addition, Dzugutov assumed that the excess entropy can be approximated by the two-body expression denoted: If the S2 description is a very good approximation for simple model systems like Lennard-Jones systems, Hoyt et al.43 have shown that it is less accurate for the density dependent many-body EAM potentials. Therefore, we test Dzugutov's idea using our computed diffusivities and radial distribution functions within the LDA approximation. Results displayed in Fig. 7 suggest that the S2 approximation is also not completely reliable when utilizing either LDA or GGA calculations.

Figure 7

The scaled diffusion coefficient vs the S2 approximation using ab initio calculations and vs the excess entropy calculated by the Carnahan Starling (CS) equation.

Thus, as for embedded atom potentials, a more reliable determination of the excess entropy is needed. The excess entropy can be expressed formally as a multiparticle correlation expansion: SE = S2 + S3 + …, where S is the entropy contribution due to n-particle spatial correlations44. Using this approach is not obvious since Yokoyama et al.45 have shown that the expansion converges very slowly for metallic systems since the four-body and even higher order terms of correlations have non-negligible contributions. Hoyt et al.43 used a lambda integration technique to get the excess entropy at T = 1500 K followed by Monte Carlo simulations to evaluate the excess entropy at any desired temperature.

Here we adopt a different strategy by making use of the Carnahan-Starling equation which is known to give the best analytical description of the thermodynamics of hard-sphere systems as long as the packing density ξ is smaller than 0.5 (Ref. 46). Note that the packing density is given by ξ = πρσ3/6. In this case, the excess entropy can be derived from the Carnahan-Starling equation as All entropy calculations can be done using the packing fraction ξ and its evolution with temperature. We determine ξ using the experimental densities8 and σ values obtained by considering that the pair-correlation entropy of the hard sphere model, , is equal to the pair-correlation entropy provided by LDA calculations, . The pair correlation function of the hard sphere model, g(r), is computed by the thermodynamically self-consistent integral equation method47. The evolution of the packing fraction with temperature is given only by the evolution of the experimental density with temperature since g(r) and then σ do not vary with temperature as obtained in our simulations. We can mention that ξ varies from ξ = 0.466 at T = 875 K to ξ = 0.415 at T = 1250 K. An important point is that these values lead to similar values of the isothermal compressibility calculated using either the Carnahan Starling expression, i.e. ρk = (1 − ξ)4/(2ξ(4 − ξ) + (1 − ξ)4), or the extrapolation of S(q) at q = 0, i.e. ρk = S(0), S(q) being obtained from LDA calculations. Table I displays the close correspondence between the two sets of calculations, which is an indication of the thermodynamic consistency of our approach to determine the excess entropy.

Table 1

values of the packing fraction ξ and isothermal compressibility χT as a function of temperature T

T (K)	ξ = πρσ3/6	χT = (1 − ξ)4/(2ρkBTζ(4 − ξ) + (1 − ξ)4)	χT = S(0)/ρkBT
875	0.466	0.024	0.029
1000	0.446	0.029	0.031
1125	0.431	0.033	0.034
1250	0.415	0.038	0.040

Using this simple analytical expression of the excess entropy, we plot in Fig. 7 the relationship between the excess entropy and the dimensionless form of the diffusion coefficient, D*. Clearly the Dzugutov scaling law becomes valid. We plot also the relationship using Eq. (10) with . The results clearly indicate a necessity of introducing the temperature dependence of ξ or σ, when the hard sphere model is employed.

In summary, we have computed structural and dynamic properties of liquid aluminum by ab initio molecular dynamics using two different exchange and correlation potentials, LDA and GGA. We determine the temperature dependence of the self-diffusion coefficient as well the viscosity using direct methods. Comparisons with experimental data favor the LDA approximation to compute dynamic properties of liquid aluminum. We show that the GGA approximation induce more important backscattering effects due to an enhancement of the icosahedral short range order (ISRO) that impact directly dynamic properties like the self-diffusion coefficient.

Our findings also support the validity of using the SE relation to describe dynamic properties of liquid Al on the basis of ab initio molecular dynamics calculations.

Finally, we show that Dzugutov's scaling law, which relates the scaled diffusivity to the excess entropy, is obeyed provided the excess entropy as computed accurately. Our results indicate that the simple two-body S2 approximation is not sufficient and that the excess entropy determined by the Carnahan Starling equation, including temperature evolution of the packing fraction, is adequate. Note that the excess entropy based on the Carnahan Starling equation is given by an analytical expression which can be easily extended to binary mixtures. Such a work is under consideration.

Methods

The AIMD simulations were performed within the framework of plane-wave-based density-functional theory as implemented in the Vienna ab initio simulation package48. A plane-wave basis set with an energy cutoff of 241 eV was used. The electron-ion interaction was described by the projected augmented-wave method49, and both the local-density approximation20 and the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) form2122 for the exchange-correlation energy were investigated. All the dynamical simulations were carried out in the NVT ensemble by means of a Nosé thermostat to control temperature. Newton's equations of motion were integrated using Verlet's algorithm in the velocity form with a time step of 1 fs. A cubic unit cell with periodic boundary conditions containing 256 atoms was used in the simulation. Only the Γ-point sampling was considered to sample the supercell Brillouin zone. The liquid samples were first prepared at 1500 K to reach thermal equilibrium (well above the melting temperature of 933 K) followed by cooling to 1250, 1125, and 1000 K successively with a rate of 1013 K/s. We also explore the undercooling state by performing runs at 875 K (LDA) and 850 K (GGA). At each temperature, the volume V of cell was chosen to reproduce the experimental densities8 and after an equilibration for 3 ps, the run was continued for 80 ps. The remainder of the simulation was used to extract the physical properties. 2000 configurations were used to produce averaged quantities such as the static structure factor and the pair-correlation function. Among these configurations, we have selected ten independent ones, regularly spaced in time (each separated by a time interval close to 8 ps), to extract their inherent structures50. To this end, the steepest-descent energy-minimization procedure with the conjugate gradient method is imposed on each of these configurations. This method allows us to uncouple the vibrational motion from the underlying structural properties, since atoms are brought to a local minimum in the potential-energy surface.

Author Contributions

N.J. and A.P. conceived the research, performed simulations discussed and wrote the manuscript. The authors declare that they have no competing financial interests. Correspondence and requests for materials should be addressed to N.J.