Ab-initio free energies of liquid metal alloys: application to the phase diagrams of Li-Na and K-Na

Comparison of free energies between different phases and different compositions underlies the prediction of alloy phase diagrams. To allow direct comparison, consistent reference points for the energies or enthalpies are required, and the entropy must be placed on an absolute scale, yielding absolute free energies. Here we derive absolute free energies of liquids from ab-initio molecular dynamics (AIMD) by combining the directly simulated enthalpies with an entropy derived from simulated densities and pair correlation functions. As an example of the power of this method we calculate the phase diagrams of two binary alkali metal alloys, Li-Na and K-Na, revealing a critical point and liquid-liquid phase separation in the former case, and a deep eutectic in the latter. Good agreement with experimental data demonstrates the power of this simple method.


I. INTRODUCTION
The free energy of a compound, G = H − T S, depends on both the enthalpy H and the entropy S. Shifting H by a constant value has no consequence because only free energy differences enter into thermodynamics. Although the values are arbitrary, comparing the free energies of competing phases is facilitated if consistently chosen reference values are applied to the enthalpies. By thermodynamic convention, the reference point for enthalpy of a compound is chosen as the sum of the enthalpies of all constituent elements in their stable state at standard temperature and pressure. However, other choices can be equally valid, in principle. For example, one could choose the enthalpy as calculated within density functional theory at temperature T = 0K as the reference point, placing high temperature enthalpies as calculated from ab-initio molecular dynamics (AIMD) on a well-defined absolute scale. In contrast, since the entropy enters the free energy multiplied by T , its actual value cannot be chosen arbitrarily. Indeed, a unique reference point for the entropy is provided by the Third Law, namely S must vanish at T = 0K. The combination of enthalpy on a well-defined relative scale and entropy on an absolute scale yields absolute free energy.
Knowledge of absolute entropy S(E) is equivalent to knowledge of the configurational density of states Ω(E). Histogram 1 and entropic sampling methods 2 such as  calculate Ω(E) up to an unknown constant factor, yielding relative but not absolute entropies. In some special discrete cases, such as lattice models where the total number of states is known, Ω(E) can be normalized yielding absolute entropy.
Most methods to compute absolute free energy rely on connecting the free energy of interest to some reference state of known free energy. Thermodynamic integration 4 calculates ∂G/∂λ through simulation (λ is some parameter in the Hamiltonian) then numerically integrates this derivative. Thermodynamic perturbation theory 5,6 expresses G(λ) as a low-order Taylor series expansion. Umbrella sampling 7 and the Bennett acceptance ratio method 8 provide increased computational efficiency to these basic approaches, as do other schemes 9 .
In practice, the thermodynamic integration requires reversible paths so that the free energy and its derivatives are well-defined. This is sometimes referred to as the slow growth approach 10-12 . Jarzynski's identity exp(−β∆F ) =< exp(−βW λ ) > holds even for nonequilibrium transitions and allows for evolution over short time durations, which is known as the "fast growth" method 13 .
2 Other methods build the free energy through sequential addition of particles. The exact scanning approach [14][15][16] computes the partial density of states ρ(α k |α k−1 · · · α 1 ) where the set {α j } is an ordered sequence of states containing successively more particles. Similarly, the particle insertion method 17 calculates the free energy difference of a k-atom system and a (k − 1)-atom system (i.e. the chemical potential), although without explicitly calculating the density of states.
The empirical CALPHAD approach 18,19 , proposes analytical free energy models for the Gibbs free energy G(x, T ) of a compound with composition x. Starting from the ideal free energy G ideal = H ideal − T S ideal , CALPHAD models the excess free energy in a series of Redlich-Kister polynomials 20 with coefficients obtained from experimental information such as heat capacity and phase diagrams. The result is a set of free energy functions in analytic form that can be used to interpolate the free energy into compositions for which no data is available.
We recently developed an approximate method to calculate absolute entropy of liquid metals from AIMD simulations of their densities and pair correlation functions 21,22 . Since entropies cannot be derived directly from simulations, our method provides a feasible approach to calculate the absolute entropy, and hence the absolute free energy, with high accuracy and reduced computational effort. As a demonstration of the utility of absolute free energy and the power of our calculational approach, we apply the method to calculate the free energies of two binary alkali metal alloy systems, Li-Na and K-Na. These examples are chosen because, despite the seeming chemical similarity, the two systems exhibit very different phase behaviors. Li and Na are nearly immiscible, both in the solid state and in the liquid below a critical point at T = 578K. In contrast, K-Na compounds remain liquid below the melting points of elemental Na and K, forming a deep eutectic at 33% Na and T = 260K (−13C). The low melting point of K 2 Na (often abbreviated as "NaK") makes the liquid alloy useful as a coolant for nuclear reactors 23 and other applications. The immiscibility of Na in Li makes it potentially useful for supression of dendrites in Li-ion batteries 24 .
In the following we first describe our simulation methods, including the manner in which we obtain absolute enthalpies and entropies. Then we validate the methods by comparing our calculated densities for pure elemental liquids with experimental values, and our absolute entropies with values tabulated in the NIST-JANAF tables 25 . Finally, we present our predicted phase diagrams of Li-Na and K-Na and compare with published experimental 3 results and they show a good agreement. We find a high positive energy of mixing between Li and Na atoms which leads to phase separation at moderately high temperatures, with a critical point for phase mixing at higher temperatures. For the K-Na system, the energy of mixing is still positive, but it is relatively weak. In consequence the entropy dominates the free energy, and after incorporating the Gibbs free energies of competing solid phases we observe a deep eutectic transition at temperatures below 0C.

II. METHODS
We simulate liquid K-Na and Li-Na using ab-initio energies and forces to accurately reproduce their configurational ensembles. Different strategies are used for elemental and binary metallic systems. For pure elements standard ab-initio molecular dynamics (AIMD) simulations are performed. For binary alloys,we supplement AIMD with additional Monte Carlo chemical species swapping steps 26 in order to accelerate the sampling of diverse configurations. Enthalpies are taken directly from the ab-initio total energies, while entropies are obtained from integrals of correlation functions 21,22 . We carry out our simulations in canonical ensembles although our entropy model is expressed in the grand canonical ensemble, relying on locality of the correlations to achieve ensemble independence 27 .
Specific simulated temperatures and compositions are chosen to cover the relevant soluble regions of the Li-Na and K-Na phase diagrams. Only soluble phases are chosen to avoid contaminated our data with multiple phases and interfacial free energies. We then fit the Gibbs free energy to an analytical model from which we derive the phase diagrams by computing the convex hull of G(x, T ). Individual data points are given in the Supplemental

Material.
A. Ab-initio molecular dynamics and Monte Carlo Our AIMD simulations apply electronic density functional theory as implemented in the Vienna Ab-initio Simulation Package (VASP 28,29 ). First-principles energies and forces are calculated using the PBE generalized gradient approximation 30,31 . MD time steps are set at 1fs with the temperature controlled in the N V T ensemble using a Nose thermostat. We take a plane-wave-basis set with a cut-off energy of E = 300eV. Semi-core electrons are included 4 in the pseudopotentials for potassium atoms and sodium atoms while only valence electrons are considered for lithium atoms. We employ simulation cells of 300 atoms for K-Na and 500 atoms for Li-Na. Justification for these decisions is presented in Appendix A.
Simulations at a given temperature and composition are pre-annealed for a minimum of 1 ps, until the onset of equilibrium energy fluctuations, followed by data acquisition for a minimum of 2 ps. We take equilibrated configurations from high temperature runs as initial conditions for lower temperature runs. To predict the density at a given temperature and composition, we monitor the total pressure at five different volumes, and then find the volume at which the pressure vanishes by fitting to a quadratic.
In binary systems, we additionally perform Metropolis Monte Carlo by testing a randomly chosen interchange of two atoms' chemical species and accepting the change with probability exp (−∆E/k B T ). We attempt one species swap every 10 MD steps. On average, a total number of 300 atomic swap attempts are made with a acceptance rates around 15%.
Supplementing AIMD with Monte Carlo (MCMD 26 ) accelerates the approach to equilibrium and enhances the configurational diversity of the simulated ensemble, as discussed in Appendix B.

B. Entropy
We calculate absolute entropies directly from MCMD simulations performed at the temperatures, densities and compositions of interest by evaluating the leading terms in an expansion of the entropy in a series of progressively higher-order correlation functions 27,[32][33][34] . This method has been previously validated for elemental liquid Al and Cu, and applied to the AlCu binary liquid alloy 21,22 .
The quantum mechanically derived absolute entropy of the ideal gas is where α denotes atomic species, λ α = h 2 /2πm α k B T is the thermal de Broglie wavelength of species α, x α is its fraction, and ρ is total atomic number density. Note that S Ideal contains the entropy of mixing  The leading term in the correlation function expansion is the single-body entropy S 1 = S Ideal − k B , with the difference arising from interchange of multiple atoms 22,27 . The twobody corrections to S 1 include a fluctuation term that is positive but very small (it is proportional to the isothermal compressibility) and an information term that is negative-definite and reflects the entropy reduction due to the information content of the pair correlation functions. We approximate the total entropy as S ≈ S 1 + S 2 with S 2 = S fluct +S info . Note that S 2 can be decomposed into partial contributions for each species reflecting the loss of entropy due to the correlations. The entropy loss is larger for Li 2 Na than for K 2 Na because of the strong chemical order that prefers like neighbors (Li-Li and Na-Na), while this effect is nearly absent in the case of K 2 Na.

C. Electronic free energy
Electronic free energies are included in every case to supplement our systematic calculations of thermodynamic quantities of alkali metals. Electronic free energies at finite temperature are obtained from the DFT-predicted electronic density of state D( ) as discussed in 37 . Electronic contributions to the relative free energies ∆G are relatively small compared with ∆G itself.

D. Interpolation
Because we carry out simulations at discrete temperatures and compositions, but we wish We fit the excess enthalpy ∆H(x, T ) and two-body term S 2 (x, T ) to a quartic polynomial The simple linear temperature dependence approximation is designed for accuracy over a narrow temperature range. This approximation works well for single species liquid enthalpies as shown in Fig. 2. The S 1 term is calculated from Eq. 1 using a quartic function to fit the composition-dependent density. By this approach, we capture the logarithmic singularities of the entropy near x = 0 and x = 1.

E. Solid phases
The Gibbs free energies of competing phases must be included to determine the global phase diagrams. These phases are body centered cubic solid phases of pure elements and a Laves phase (Pearson hP12, Strukturbericht C14) KNa 2 binary phase. The Gibbs free energy for a solid phase includes the vibrational free energy G v , the electronic free energy G e , and a configurational free energy G c . In principle the configurational term includes contributions due to chemical and vacancy disorder 37 , however experimental evidence suggests that the K-Li-Na solid phases are nearly stoichiometric, so we simply approximate G c with the enthalpy of the fully relaxed (i.e. T =0K) structure. First-principle vibrational free energy calculations use the same pseudopotentials and exchange-correlation function as for the liquid simulations, but with an increased plane wave cut-off energy of 500eV and an increased electronic k-mesh density so that we may obtain accurate interatomic force constants. The Gibbs free energies at finite temperatures are calculated within the quasi-harmonic approx-  Positions are plotted from back to front, with diameters indicating depth.
the electronic entropy provides an excellent match to experiment at temperatures above melting. A small deficit remains in the liquid entropy close to the melting point that is presumably due to three-and four-body correlations 22 .

B. Binary K-Na
As shown in our simulation snapshot of K 2 Na at T=473K in Fig. 3, K and Na atoms are uniformly distributed, indicating a homogeneous liquid state. We confirm a stable mixture of liquid K-Na alloy at T=473K and x = 33%. The stability of the liquid alloy is confirmed by our calculated enthalpies, entropies and Gibbs free energies as plotted in Fig. 4. The internal energy of compound is positive, which could indicate phase separation, but the amplitude of its peak is substantially below the amplitude of the peak in the Li-Na system, and far below the magnitude of T S. Values of S 1 contain the density ρ, which varies monotonically and smoothly with composition, and also contains the entropy of mixing S Mix . Hence −T ∆S 1 is negative and strongly convex. Note that although S 2 is negative-definite, we find that  −T ∆S 2 , defined relative to the pure elements, is negative and seemingly is also convex.
As a result, the total Gibbs free energy is dominated by entropy and is convex over all compositions, resulting at continuous miscibility of K and Na at 473K. We would expect phase separation below T=189K based on extrapolation of ∆G to low temperatures, but this is preempted by the eutectic transition to the solid phases, as we now discuss.
We compare the free energy of the liquid phase with its competing solid phases in order to predict the composition-and temperature-dependent K 1−x -Na x phase diagram as shown in Fig. 4 (d). K-Na has three known low-temperature phases-elemental K (x = 0), elemental Na (x = 1) and the KNa 2 Laves phase (x = 2/3). In the experimental phase diagram, a deep eutectic transition occurs at x E ≈ 33% and T E ≈ 260k, where the K-Na alloy exists in the liquid state at temperatures below the melting points of elemental K and Na. In our predicted phase diagram, a eutectic transition is found near T E = 268K and x E = 40.3%, not very far from the experimental position.
The deviation of these two transition point might be an effect of the systematic error in prediction via DFT or it might be due to our approximations for the entropy. To understand which is most responsible, we compare the predicted melting points (i.e. the temperatures at which solid and liquid free energies cross) with experiment as seen in Fig. 4

C. Li-Na
We present a snapshot of a 300-atom Li 2 Na system at T=473K in Fig. 3(b). The snapshot shows clear phase separation into an Na-rich region and an Li-rich region that mutually coexist in equilibrium. The separation is also evident in the correlation functions in Fig. II B(b), where the amplitude of the first two peaks of the Li-Na pair correlation function is substantially smaller than those of Li-Li and Na-Na, indicating effective Li-Na repulsion.
In order to understand how internal energies and entropies contribute to the total free energies and drive the system toward phase separation, we plot our calculated enthalpy, entropy and Gibbs free energy for Li 1−x Na x in Fig. 4b. The energy cost of mixing Li and Na is large compared to that of mixing K and Na (Fig. 4a), and similar in magnitude but opposite in sign to −T S. The resulting free energy lacks convexity and hence explains the separation of liquid Li and Na at low and moderately high temperatures. Specifically, G(x) lies above its own convex hull over the interval from x = 0.07 to x = 0.90 (see dotted orange line in Fig. 4b). A liquid alloy in this composition range will phase separate into a mixture of those two endpoint compositions. Note that we only use data at compositions that lie within our predicted single phase regions.
Collecting data similar to that of Fig. 4b at higher temperatures (see Fig. C in Appendix C), we then fit the temperature evolution of the phase boundaries. Fig. 4d  In addition to the Li-Na phase separation, a eutectic transition (not shown) occurs at x Na ≈ 97% and T ≈ 290K and a monotectic transition near the melting point of Li. Both features are also reported in the experimental phase diagram 45 .

IV. CONCLUSIONS
In summary, we systematically studied thermodynamic properties of solid and liquid K-Na and Li-Na metallic alloys at finite temperature and zero pressure. The Gibbs free energies of BCC K, Li and Na, and Laves phase KNa 2 , were calculated in the quasi-harmonic approximation as implemented in Phonopy 38,39 . Standard ab-initio molecular dynamics and Monte Carlo/molecular dynamics simulations modeled the liquid alloys. Absolute entropies in the liquid state were obtained as functionals of simulated densities and pair correlation functions. We note that a similar approach is possible in the solid state also 21,46 . At T =473K we observed phase separation in Li 2 Na in contrast to phase mixing in K 2 Na, and these observations were justified by explicit calculation of the composition-dependent absolute Gibbs free energy G(x) that revealed nonconvexity in the case of Li-Na.
Extending our calculations to other temperatures, we predicted composition-temperature phase diagrams that agreed well with experiment in most respects. Specifically, we obtained a deep eutectic transition in K-Na and liquid-liquid phase separation in Li-Na. Our principal shortcoming was our inability to accurately model the critical point for Li-Na phase our results, seeking to achieve a balance between accuracy and computational efficiency.
Because the enthalpy and entropy both depend sensitively on the density of the liquid, and experimental data on density is readily available and presumably reliable, we take the accuracy of our predicted density as a measure of accuracy overall. When we simulated the Li-Na and K-Na binaries, we found that the N = 300 atom cells

Cutoff energy
The incompleteness of the plane wave basis set creates systematic errors in the calculated pressure (the Pulay stress 48 ) that diminish as the plane wave energy cutoff increases. Different elements, and even different pseudopotentials for the same element, have very different default energy cutoffs. The defaults are 140.000 eV for "Li 17Jan2003", 116.731 eV for "K pv 17Jan2003" and 259.561 eV for "Na pv 19Sep2006". The " pv" subscripts indicate that p semicore electrons are treated as valence. Calculational cost grows as the cube of the number of plane waves, which itself grows as the 3/2 power of the energy cutoff, leading to rapid growth of cost vs. cutoff, yet consistency of calculated energies requires that the applied energy cutoff be uniform across different compositions. Since the Na potential requires the highest energy cutoff for both the Li-Na and K-Na alloy systems, we explore the sensitivity of the density of Na to the energy cutoff. Testing values 260, 300 and 340 eV, with systems of N = 300 atoms at T=473K, we obtained densities of 0.02948, 0.02378 and 0.02381 atoms/Å 3 . Thus we settle on a cutoff of 300 eV. Note the experimental value is around 0.0237 (see Fig. A 2)  maining slightly high.
Note that the predicted density of K drops with increasing energy cutoff, and converges to nearly 5% below the experimental value at this cutoff. The situation is marginally improved if we switch to the "K sv" potential, but not sufficiently to compensate for the increased electron count.

XC functional
The choice of exchange correlation functional can lead to systematic errors in the density.
We compared the local density approximation (LDA) with the PBE generalized gradient approximation and with PBEsol for a system of N = 300 Li atoms at T=473K using an energy cutoff of 300 eV and a single k-point. These yielded densities of 0.04535, 0.04516 and 0.04511 atoms/Å 3 , for LDA, PBE and PBEsol respectively, compared with experimental values that range from 0.0444 to 0.0448. All three overestimate the density. PBEsol proves only marginally better than PBE, and we prefer PBE because it is more widely used. To test if this is caused by neglect of the Li core electrons, we tested the valence-3 "Li sv 10Sep2004" potential at energy cutoff 650 eV (30% above its default of 499.034 eV) and found density 0.0462, which is far above the experimental value. Thus we settle on the PBE functional and stick with our decision to use the valence-1 Li potential.

Uncertainties in thermodynamic quantities
Our calculations are subject to both systematic and statistical errors. Density functional theory itself relies on the approximate exchange-correlation functional, leading to systematic errors that we do not attempt to quantify, beyond noting that differing choices of functional had modest influence on the density as discussed in Sec. A 3, while DFT may lead to errors in the vicinity of 15K for the melting point of Na as discussed in Sec. III B. We truncate our expansion of the entropy at the pair level, systematically omitting three-and higherbody correlations. This may be the reason that our entropy falls below experiment close to the melting point 22,46 , as seen in Fig. 2. We restrict our simulations to certain finite sizes, leading to errors on the density as shown in Fig. A 1, but also leading to premature phase separation in Li-Na even at high temperatures above the critical point, when the correlation length grows beyond our simulated cell size (see Fig. C).
To estimate our statistical errors, we break our runs into three segments and evaluate the standard error on the assumption of uncorrelated errors. For K 2 Na at T=473K we find statistical errors of order 1 meV/atom in both the enthalpy H and the entropy −T S 2 .
However these are anticorrelated so that the statistical error on G is of order 0.4 meV/atom.

SUPPLEMENTAL DATA
The following are individual calculated data points for Fig. 4 and Fig. C