Yiheng
Shen
and
Wei
Xie
*
Materials Genome Institute, Shanghai University, Shanghai 200444, China. E-mail: xiewei@xielab.org
First published on 18th February 2025
Understanding the phase stability of elemental lithium (Li) is crucial for optimizing its performance in lithium-metal battery anodes, yet this seemingly simple metal exhibits complex polymorphism that requires proper accounting for quantum and anharmonic effects to capture the subtleties in its flat energy landscape. Here we address this challenge by developing an accurate graph neural network-based machine learning force field and performing efficient self-consistent phonon calculations for bcc-, fcc-, and 9R-Li under near-ambient conditions, incorporating quantum, phonon renormalization and thermal expansion effects. Our results reveal the important role of anharmonicity in determining Li's thermodynamic properties. The free energy differences between these phases, particularly fcc- and 9R-Li are found to be only a few meV per atom, explaining the experimental challenges in obtaining phase-pure samples and suggesting a propensity for stacking faults and related defect formation. fcc-Li is confirmed as the ground state at zero temperature and pressure, and the predicted bcc-fcc phase boundary qualitatively matches experimental phase transition lines, despite overestimation of the transition temperature and pressure slope. These findings provide crucial insights into Li's complex polymorphism and establish an effective computational approach for large-scale atomistic simulations of Li in more realistic settings for practical energy storage applications.
Previous first-principles density functional theory (DFT) simulations have explored the energetics of various polymorphs of Li, predicting near degeneracy between bcc- and fcc-Li.20 The inclusion of zero-point energy and finite-temperature anharmonic effects further influences the relative stability, e.g. the vibrational entropy drives the transition from close-packed structures to bcc-Li upon heating,21 and nuclear quantum effects significantly influence the phonon spectra and phase stability.16 However, these accurate DFT-based calculations are computationally demanding, especially for large-scale simulations required for studying finite-temperature properties. To our knowledge, such limitations have led to the exclusion of the 9R-Li from accurate thermodynamic calculations that properly account for quantum and anharmonic effects, due to its complex crystal structure with significantly lower symmetry than the two competing phases of bcc- and fcc-Li.
The development of accurate and efficient interatomic potentials has been instrumental in enabling large-scale simulations of Li at length and time scales inaccessible to DFT.22 Recently, machine learning force fields (MLFF),23–25 in particular those based on equivariant graph neural networks (GNN),26–28 have emerged as a promising tool for achieving DFT accuracy in interatomic potentials, while retaining computational efficiency,29,30 offering the opportunity to investigate thermodynamics of Li with greater accuracy and computational efficiency. Wang et al.29 developed a deep learning potential (DP) for elemental Li and used it to calculate various bulk, defect and surface properties. Phuthi et al.30 developed a GNN-based MLFF in the NequIP architecture for elemental Li that accurately predicted the small difference (about 2 meV per atom from DFT) in the zero-temperature total energies of fcc- and bcc-Li, and examined elastic properties, adsorption energies and surface diffusion. However, a systematic thermodynamic analysis of Li leveraging the full power of advanced MLFF remains to be conducted.
In this work, we conduct a comprehensive study of the phase stability of Li in its bcc, fcc, and 9R phases under near-ambient conditions. We develop a state-of-the-art equivariant GNN-based MLFF in the MACE architecture27,28 trained on DFT data, and perform self-consistent phonon (SCP) calculations based on it, accounting for nuclear quantum and the anharmonic vibrational effects of both phonon renormalization and thermal expansion. Temperature and pressure dependent Gibbs free energies of the three phases are calculated afterwards based on the effective force constants obtained from SCP, facilitating the investigation of the stability of the three phases across a range of temperatures and pressures with high fidelity. Our results provide insights into the phase stability of Li and offer practical advances in modeling the thermodynamic properties of this vital element.
Using the MACE force field as the energy and force calculator, we then perform SCP calculations for bcc-, fcc- and 9R-Li at various temperatures and lattice parameters. Lattice parameters are measured in this communication by strain as calculated with respect to their relaxed values at zero temperature. For 9R-Li which has hexagonal lattice system with two independent lattice parameters a and c, a mesh of biaxial strains εa and uniaxial strains εc are considered to account for anisotropic thermal expansion. Details of the SCP calculations are summarized in Note S2 in the ESI.† In essence, the nuclear quantum effects are included by following Bose–Einstein statistics in the sampling of the configuration spaces, phonon renormalization is treated in the SCP calculations each performed at fixed lattice parameters, while thermal expansion is treated by minimizing the Helmholtz free energies from SCP calculations at different lattice parameters following the condition of equilibrium for NPT ensemble, similar to how it is treated in quasi-harmonic approximation. The MACE force field easily accelerated the above calculations by over three orders of magnitude with respect to DFT-VASP calculation, as measured by the wall time of a static calculation of one SCP supercell. See Note S2 of the ESI† for more detailed comparison of the computational cost of MACE and DFT-VASP calculations.
Fig. 2 shows the temperature dependent phonon band structures of bcc-, fcc- and 9R-Li at select lattice parameters, as measured by strain. Fig. S6–S8† provide the results at all calculated lattice parameters. The renormalized phonon modes slightly harden (i.e. the mode frequencies increase) upon heating and profoundly harden upon compression. The dependence of frequency on lattice parameters, which leads to thermal expansion, is found to be more significant than phonon renormalization due to temperature. The two anharmonic effects seem to couple more strongly when Li is subject to compression, and at the strain value of −0.2, soft modes develop for both bcc-Li and fcc-Li and significant phonon renormalization is found for the two phases. On the other hand, 9R-Li seems to be less sensitive to compression, and remains dynamically stable in the whole range of strains explored in this study. The different reactions to strain are in agreement with the fact that 9R-Li has only been observed under low pressure,18 as it gains less entropy due to anharmonicity upon increasing pressure than the competing fcc-Li phase.
An important question is then the strength of Li's anharmonicity. To provide a quantitative estimate, we calculate the anharmonicity measure43 for the three studied phases. As shown in Fig. 3, the anharmonicity measures of the three phases are comparable to one another. Those of bcc- and fcc-Li reach 0.3 at 300 K, while that of 9R-Li should still be about the same, as judged from the trend. This result indicates that about 30% of the interatomic forces of elemental Li originate from anharmonic interactions. As a reference, at 300 K the anharmonicity measure of crystalline silicon, a typical material with weak anharmonicity is about 0.2, while that of a strongly anharmonic halide perovskite KCaF3 is about 0.5.43 This result suggests that anharmonicity in Li is not negligible, which is in accordance with what can be expected from the weak metallic bonding and the light atomic weight of Li, as discussed before,16,44 and therefore warrants the careful treatment of both phonon renormalization and thermal expansion in this study.
![]() | ||
Fig. 3 Anharmonicity measure of bcc-, fcc- and 9R-Li as a function of temperature. Analogous results excluding nuclear quantum effects are provided in Fig. S11 in the ESI.† |
After the lattice parameter and temperature dependent effective force constants are obtained by SCP calculations, we can now calculate the thermodynamic functions of the three phases of Li. The Helmholtz free energy F of a phase with lattice A at a given temperature T, i.e. F(T, A), is obtained as the sum of the vibrational term Fvib calculated based on the effective force constants45,46 and the electronic term Felec calculated based on Fermi–Dirac smearing in DFT by VASP. Details of the Fvib calculations are explained in Note S2 in the ESI.† The Gibbs free energy G at given temperature T and pressure p, i.e. G(T, p), is then calculated by minimizing the weighted 3rd-order polynomial fit of F(T, A) + p·det(A) over A. We stress that this last step is a non-trivial task, and significant caution needs to be exerted to avoid errors caused by poor fitting, as explained in Note S3 in the ESI.†
The predicted temperature- and pressure-dependent Gibbs free energy and equilibrium lattice parameters for bcc-, fcc- and 9R-Li are summarized in Fig. S9 and S10,† respectively. To provide a glimpse into these full results, Fig. 4a shows that the predicted primitive-cell volume of bcc-Li coincides with the experimental values at room temperature very well, despite slight underestimation at zero-pressure, which validates the accuracy of our thermodynamic calculations. Fig. 4b shows the pressure-dependent Gibbs free energy at 200 K, in which fcc-Li is the ground state throughout the explored range of pressure in agreement with the recent experimental study.18 It is also worth mentioning that, the Gibbs free energy differences between bcc- and fcc-Li fall in the few meV per atom range at zero pressure. We recognize that such energy differences are also close to the MAE of the MACE force field, and suspect that the predicted phase diagram has inherited some inaccuracies, as will be discussed below.
![]() | ||
Fig. 4 Pressure dependent (a) equilibrium primitive-cell volume, (b) Gibbs free energy (reference to fcc-Li), (c) phase diagram and (d) Gibbs free energy difference between fcc- and 9R-Li. Experimental results in (a) are taken from ref. 47. White lines in (c) and (d) represent experimental transition lines of 7Li, showing that upon cooling at low pressure bcc-Li actually transits to metastable 9R-Li instead of fcc-Li.18 |
Finally, the continuous p–T phase diagram of Li is predicted in Fig. 4c. bcc-Li is identified as the equilibrium phase at ambient temperature and pressure and its phase boundary with fcc-Li shows an upward concave curvature with an extreme point near 0.5 GPa and 300 K. The elevated phase boundary in bcc-Li with respect to pressure above 0.5 GPa qualitatively reproduces the trend of the experimental transition lines between bcc- and fcc-Li, despite a systematic overestimation of the slope and the critical temperature value (∼300 K vs. ∼100 K from experiment18). The transition lines are only rough estimation of the phase boundary, as they could be significantly affected by various kinetic factors and also depend on the actual experimental paths taken, as discussed extensively in ref. 18. The seemly large difference between our computationally predicted phase boundaries and the experimental transition lines is also not very surprising given the very flat energy landscape of Li. In fact, Fig. 4d shows that the Gibbs free energy difference between the two close-packed stacking variants fcc- and 9R-Li in the whole 0–4 GPa and 0–250 K range is only 2–3.6 meV per atom. This means that small deviation in Gibbs free energy can lead to the prediction of significantly different phase diagram, as discussed in Note S3 of ESI.†
The findings of this study hold significant implications for advancing lithium metal battery technologies, particularly in addressing the persistent challenge of dendrite formation. The near-degeneracy in Gibbs free energy between bcc-, fcc- and 9R-Li phases in wide ranges of temperature and pressure suggests that even minor perturbations during battery operation—such as localized stress, temperature fluctuations, or interfacial interactions—could trigger phase transitions or coexistence. Such metastability may promote the formation of stacking faults, twin boundaries and other defects in bulk Li, which act as nucleation sites for heterogeneous lithium deposition and dendrite growth. This underscores the critical need to engineer electrode architectures or electrolyte interfaces that thermodynamically stabilize a single phase under operational conditions. For instance, substrate materials or artificial solid-electrolyte interphases (SEI) designed to epitaxially template bcc-Li could suppress defect formation and encourage uniform plating.48,49 Furthermore, the pronounced anharmonicity and nuclear quantum effects revealed in this work highlight the necessity of incorporating these factors into computational models to accurately predict lithium's behavior in real-world battery environments. The demonstrated success of the MACE-based approach in capturing subtle phase equilibria at a small computational cost provides a powerful tool for screening materials and operating conditions (e.g., temperature and strain) that favor phase homogeneity through large-scale atomistic simulations of Li in more realistic settings. By integrating such insights, strategies such as strain-modulated cell designs or thermal management systems could be devised to stabilize desired Li phases, mitigate stress concentrations, and ultimately enhance the cycling stability and safety of lithium metal batteries. This study thus bridges fundamental understanding of Li's complex phase behavior to actionable strategies for realizing practical high-energy-density lithium metal anodes.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4ta08860c |
This journal is © The Royal Society of Chemistry 2025 |