Prediction of stable Li-Sn compounds: boosting ab initio searches with neural network potentials

The Li-Sn binary system has been the focus of extensive research because it features Li-rich alloys with potential applications as battery anodes. Our present re-examination of the binary system with a combination of machine learning and ab initio methods has allowed us to screen a vast configuration space and uncover a number of overlooked thermodynamically stable alloys. At ambient pressure, our evolutionary searches identified a new stable Li$_3$Sn phase with a large BCC-based hR48 structure and a possible high-T LiSn$_4$ ground state. By building a simple model for the observed and predicted Li-Sn BCC alloys we constructed an even larger viable hR75 structure at an exotic 19:6 stoichiometry. At 20 GPa, new 11:2, 5:1, and 9:2 phases found with our global searches destabilize previously proposed phases with high Li content. The findings showcase the appreciable promise machine learning interatomic potentials hold for accelerating ab initio prediction of complex materials.


I. INTRODUCTION
Accurate density functional theory (DFT) approximations introduced over the past three decades have proven to be powerful practical solutions for quantum mechanical description of materials properties [1][2][3]. Reliable prediction of structures' thermodynamic stability [4][5][6] has had a particularly transformative impact on the materials discovery process, as ab initio ground state searches have been used to screen large chemical spaces and identify new synthesizable materials [4,7]. Despite the growing number of confirmed predictions [8][9][10], the pace of DFT-guided discovery is ultimately limited by the high cost of ab initio calculations.
MLP-assisted prediction of ground state structures differs notably from these materials modeling applications. While ns-long MD trajectories or a-few-nm-sized configurations cannot be checked with DFT directly, pools of candidate structures with lowest Gibbs free energy found in MLP-based searches can and should be examined with the reference method. The benefit of using surrogate models should then be measured against the net cost of their parametrization and the following DFT analysis of viable candidates. MLPs have proven to be essential for finding stable nanoparticle configurations because direct global optimizations with the DFT become impractical for cluster sizes above a few dozen atoms [31][32][33], while empirical potentials generally lack the accuracy needed to effectively guide ab initio searches. In fact, our recent large-scale comparative study showed a significant advantage of a typical neural network (NN) interatomic model over traditional potentials and led to the revision of several ab initio ground states for Au clusters with 30-80 atoms [34].
In this study, we aim to illustrate the capabilities of our recently developed framework for MLP-acccelerated structure prediction [9,46] by performing a largescale exploration of Li-Sn alloys. The motivation for (re)investigating this common binary system is twofold. First, the existence of several Li-rich compounds with high electrical conductivity and better ductility compared to other group-XIV elements make them appealing candidate materials for Li-ion battery anodes [47][48][49]. Second, the binary system has been recently investigated with the state-of-the-art ab initio prediction methods in two systematic studies and shown to host new thermo-dynamically stable compounds at ambient and elevated pressures [47,50]. The development and application of an accurate NN model for Li-Sn has allowed us to probe over a million structures across the full composition range at different (P ,T ) conditions and identify several new ground states. Particularly surprising are thermodynamically stable BCC-based Li 3 Sn and Li 19 Sn 6 phases with large prototypes that appear not to have been observed in any other binary alloys. These findings illustrate the appeal of employing MLPs to accelerate ab initio prediction of complex materials.

Neural network interatomic potential
A NN model of the Behler-Parrinello type was constructed using an automated iterative scheme implemented in our MAISE-NET framework [9] that features an evolutionary sampling algorithm to generate representative reference structures and a stratified training protocol to build multicomponent models on top of elemental ones [46]. We relied on our previously developed Li and Sn NN models [9] with 51-10-10-1 architectures and σ Li E = 2.3 meV/atom and σ Sn E = 8.3 meV/atom rootmean-square errors (RMSEs). The Li-Sn NN with a 145-10-10-1 architecture and 1,880 adjustable parameters 1 was fitted to 6,046 energy and 46,410 force data in binary structures with up to 32 atoms generated in four MAISE-NET cycles. The optimization of the binary model for 120,000 steps resulted in σ Li-Sn E = 10.2 meV/atom and σ Li-Sn F = 49.0 meV/ A RMSEs over a test set of 671 structures. All elemental and binary models were based on our stardard set of Behler-Parrinello symmetry functions with a 7.5 A cutoff radius [34]. In total, the extension of the elemental Li and Sn models to the binary one cost under 35,000 CPU hours.

Structure prediction strategy
Overview Identification of ab initio Li-Sn ground states at different (P, T ) conditions included evolutionary searches with the NN model at T = 0 K, re-examination of viable T = 0 K candidates with the DFT, and investigation of phase stability at elevated temperatures at the NN and DFT levels. The global search for stable Li-Sn alloys at T = 0 K and P = 0 GPa or P = 20 GPa was driven by an evolutionary algorithm implemented in MAISE [9]. Detailed in our previous studies [9,63,64], this global optimization engine has been used primarily in combination with DFT to predict or solve complex ground states. It has led to the discovery of several confirmed new prototypes with large unit cells: oP10-FeB 4 [63,65], mP20-MnB 4 [4], and tI56-CaB 6 [64] found without any structural input as well as mP24-Cu 2 IrO 3 and tP56-Na 3 Ir 3 O 8 found by seeding the searches with known structures [66,67]. Noteworthy features of the present NN-based evolutionary searches and post-search analysis are described below.
Selection of compositions Exhaustive screening of relevant stoichiometries is a requisite for identification of stable compounds. It has been demonstrated [68,69] that unsupervised variable-composition evolutionary searches can result in an efficient location of stable stoichiometries. One of the recent studies on Li-Sn [47] employed such an algorithm and refined solutions with follow-up evolutionary searches at select compositions. Given the richness of the Li-Sn phase diagram and the complexity of the known ground states at non-trivial compositions, we relied on a supervised selection of stoichiometries for our fixed-composition searches. Figure 1 summarizes available information on Li-Sn and related chemical systems helpful for determining compositions pertinent to this study. The three lowest sets list previously synthesized Li-Sn, M -Sn (M = Na-Cs), and Li-X (X = Si, Ge, and Pb) alloys. Among the M -Sn binaries, Li-Sn stands out as a system with multiple Li-rich compounds and unique stable compositions across the full range. It is evident that the size of the alkali metals is a dominant factor defining stable M :Sn ratios, which renders the observed Na-Cs tin alloy compositions not particularly relevant for Li-Sn under ambient conditions. The Li-X set, on the other hand, provides important clues regarding possible compositions and morphologies that could be observed in the Li-Sn binary. As discussed in the following sections, some of the previously predicted compounds displayed near the top  [62]. The two grey groups specify stoichiometries observed in related M -Sn and Li-X binaries [62]. The blue set refers to thermodynamically stable compounds predicted in previous studies [47,50] for 0 GPa and 20 GPa (the latter are marked with asterisks). The top red group corresponds to additional compositions considered in this study: all possible m : n ratios with (m + n) ≤ 8 and a few larger subsets with m = 6-23 and n = 1-9 near regions with known stable compounds. of Fig. 1 were shown to be stable in these prototypes.
Our main focus was on compositions with previously reported analogs in Fig. 1 containing fewer than 40 atoms in the primitive unit cells. The very large observed cF420-Li 17 Sn 4 prototype was examined directly with the DFT, as finding global minima from scratch for crystalline structures above about 40 atoms is challenging even with evolutionary searches based on surrogate models. We also considered all possible m : n ratios with m + n ≤ 8 and included a few larger m : n subsets with m = 6-23 and n = 1-9 to sample the Li-rich end of the phase diagram (the top row in Fig. 1).
Evolutionary search settings At each selected composition, we ran separate evolutionary runs with different numbers of formula units. Typically, structures had between 1 and 8 formula units and did not exceed 40 atoms in the primitive unit cell. Randomly generated populations of 32 members were evolved for up to 300 generations with standard evolutionary operations. Namely, 8 new members were constructed from a single parent structure through mutation (random atom displacements, atom swaps, and unit cell distortions) while 24 offspring were created from two parents through crossover (combination of two roughly equal parts obtained with planar cuts) [63]. Child structures were locally relaxed with the NN potentials for up to 300 Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimization steps and assigned a fitness based on the final enthalpy. Our fingerprint method based on the radial distribution function (RDF) [9,64] was used to identify similar structures and decrease their survival probability.
Post-search DFT analysis at T = 0 K Upon completion of a NN-based evolutionary run, a subset of viable candidates from all generated local minima was selected for further re-examination with the DFT. The global DFT minimum is captured provided that it is (i) at least a local minimum on the NN PES; (ii) visited by the search algorithm; and (iii) included in the analysis pool [34]. The main adjustable parameter in the selection process is the enthalpy window. By default, we collected all distinct minima within 20 meV/atom (≈ 2σ Li-Sn E ) above the NN ground state to account for the NN typical errors in the evaluation of relative enthalpies but reduced the window at some compositions with a large number of competing low-symmetry states. We relied on a 0.95 RDF-based dot product cutoff to exclude similar structures and used a 0.1 tolerance to symmetrize the unit cells [9] before relaxing them with DFT.
Analysis of phase stability at high T Identification of high-T ground states is commonly done by including the vibrational entropy contribution to the Gibbs free energy for a pool of low-enthalpy phases found in global structure searches at T = 0 K [2,9]. The relative change in Gibbs free energy can reach a few dozen meV/atom for structures with substantially different phonon density of states (DOS) observed, e.g., in α-Sn and β-Sn [70][71][72]. Given the demanding nature of phonon calculations at the DFT level, we used the following protocol to examine high-T phase stability. Firstly, we calculated thermal corrections with the NN model (∆F NN vib ) for all distinct candidate structures in the T = 0 K DFT pools within a 20 meV/atom window. Secondly, we constructed a new convex hull at 600 K using G DFT+NN = H DFT + ∆F NN vib and selected up to 5 structures per pool with Gibbs free energies no further than 10 meV/atom from the tie-line. Finally, we evaluated thermal corrections for these structures with DFT and determined high-T ground states by examining convex hulls at all temperatures in the 0-800 K  An overview of the NN-guided identification of zero-T and high-T ab initio ground states. The flowchart illustrates the selection of viable candidates following a representative evolutionary search performed for Li18Sn4 unit cells at 0 GPa. The two sets on the left correspond to the lowestenthalpy NN minima (in red) and the resulting reoptimized DFT minima (in black). The two subsets on the right show the Gibbs free energies in the 0-600 K range with the phonon entropy contributions calculated with the NN model (in red) and the DFT (in black). The table at the bottom shows the tally of the structures examined and computational resources used in the full exploration of the Li-Sn binary.
range, in intervals of 10 K. Both NN and DFT phonon calculations were performed with the finite displacement method as implemented in Phonopy [73]. We used symmetry-preserving expansions of primitive or conventional unit cells to generate supercells with 72-216 atoms. Structures included in the thermodynamic stability analysis were dynamically stable, and any structures which remained dynamically unstable were disregarded if a new structure could not be created from their imaginary frequency eigenmodes.

A. Overview of global structure search results
The information presented in Figs. 2-4, serves to illustrate the size of the configuration space sampled in this work and the reliability of the NN-based structure prediction approach employed for identification of stable Li-Sn phases.
According to our results, populations in the evolutionary searches typically stopped improving after 210, 910, and 2,300 configurations were probed in runs with up to 14 atoms, 15-24 atoms, and over 25 atoms per unit cell, respectively. While the stochastic searches are not guaranteed to converge to the global minimum in every in-stance, we find that sampling 4,400 candidates structures with up to 24 atoms is generally sufficient. Overall, more than 1.1 million local optimizations with the NN model were performed for binary structures with no initial symmetry and an average of 19.4 atoms per unit cell, which required 952,000 CPU hours (see Fig. 2). An equivalent set of DFT calculations would require an estimated 100 million CPU hours (the increase in the average computational cost was estimated by relaxing 20 structures with 5-30 atoms for the same number of steps in the NN and DFT optimizations). The selection and symmetrization of configurations with low NN enthalpies narrowed the pool down to 4,450 structures (about 13 structures per pool) and cost considerably less to re-optimize with DFT. Importantly, the near-stable structures were described more accurately with the NN model (6.9 meV/atom RMSE compared to σ Li-Sn E = 10.2 meV/atom) and underwent relatively small geometrical changes in the following DFT local relaxations (2.0 meV/atom enthalpy decrease on average). As a result, the best DFT minima ranked among the most favored ones within the NN pools (the structures with the lowest NN enthalpy remained favored at the DFT level in 45% of all considered cases) and had either matching or lower enthalpies compared to all previously reported Li-Sn phases at T = 0 K.
Our analysis of phase stability at elevated temperatures benefited similarly from the use of the NN model. Following the selection scheme outlined in Section II, we performed phonon calculations for 2,300 and 233 structures at the NN and DFT levels, respectively. With the free energy contributions averaging 157 meV/atom at 600 K, the correspondence between the two methods was 7.2 meV/atom. With the RMSE of 1.5 meV/atom and the 4.9 meV/atom average change in relative stability at T = 600 K (see Fig. S1 for typical results), the NN-based screening allowed us to examine an order of magnitude more candidate structures and uncover possible high-T ground states at 0 GPa and 20 GPa.
It is important to note that while DFT systematic errors in calculations of thermodynamic stability may reach a few dozen meV/atom in certain cases [63,74], relative energies evaluated for similar structures tend to be far more accurate due to cancellation of errors [36,75]. Our tests performed for different functionals and convergence settings (Tables S1 and S2) demonstrate that the relative enthalpies can be resolved to within a fraction of 1 meV/atom for most considered structures with underlying BCC morphologies. Temperature estimates for transitions between near-degenerate phases with small differences in free energy corrections can have significant uncertainties of a few hundred Kelvin [36,76].
Our main findings regarding competing Li-Sn phases at 0 GPa and 20 GPa are summarized in Figs. 3 and 4. The DFT results allowed us to refine known T = 0 K convex hulls (top), identify additional metastable alloys (middle), and establish temperature ranges of phase stability (bottom). Namely, two new ambient-pressure ground states are added next to the previously pro- posed ones [47,50] and several new high-pressure ground states substantially redefine the Li-rich end of the convex hull predicted in Ref. 47. The collection of local minima obtained in our NN-based evolutionary searches illustrates the presence of many more phases within the 20 meV/atom window than found with the DFT-based random sampling approach [50]. The inclusion of thermal corrections indicates that only a few considered metastable phases have a chance of becoming thermodynamically stable.

B. Stability and morphology of Li-Sn phases
The variety of observed Sn alloy morphologies [77][78][79][80] can be attributed to the particular position of the element in the periodic table. Similarly to the C-Ge group-XIV members, Sn displays a propensity to forming covalent bonds and crystallizes in the open diamond structure (α-Sn) at low (P, T ). In contrast to the light isoelectronic elements, the more pronounced s-p relativistic splitting reduces the favorability of the s-p hybrid orbitals and, consequently, Sn adopts various ground states with metallic bonding, such as β-Sn above 286 K [70] or closed BCT structure between 11 GPa and 32 GPa [77,79]. Li is an archetypical sp metal with the BCC ground state under ambient conditions transforming into FCC above ∼ 5 GPa [81,82]. In combination, Li and Sn have been observed to form BCC-based alloys at high Li concentrations and compounds with covalent frameworks at high Sn concentrations. The prevalence of BCC motifs can be viewed as a consequence of the negative chemical pressure induced by the larger Sn (a Sn BCC = 3.81 A) that further favors the BCC-Li lattice (a Li BCC = 3.43 A) over the high-pressure FCC-Li polymorph. The largest magnitude of the formation enthalpy achieved around the 4:1 composition (Figs. 3 and 4) is consistent with the octet rule [47,[84][85][86].
At ambient pressure, Li 17 Sn 4 is the most Li-rich synthesized alloy in this binary system. Goward et al. performed a comprehensive XRD analysis of Li 22 X 5 (X = Ge, Sn, and Pb) compounds and determined that the slight compositional variability across the series originates from selected occupation of Li sites [87]. The ordered cF420 (F43m) structural model used to simulate Li 17 Sn 4 [47,50] has only one Sn Wyckoff site (16e) with the full (8+6)-atom BCC coordination that produces tetrahedral agglomerates of Sn polyhedra (see Fig.  5 (a)). The other three Sn Wyckoff sites, 16e, 24g, and 24f , are surrounded by 13 Li atoms and generate isolated Sn atoms sufficient to fill up all space with interlocking polyhedra. As in previous studies [47,50], the phase is found to be thermodynamically stable in our DFT calculations.
The observed oS36-Li 7 Sn 2 phase represents a regular BCC lattice with fully occupied sites. The Sn atoms in this Li 7 Ge 2 prototype are not fully segregated, making one of the two Sn Wyckoff sites have a Sn neighbor in , and predicted in this study (g,l). The Li and Sn atoms are shown in green and grey, respectively, while the shaded polyhedra, sections, and planes highlight distinctive morphological features discussed in the text. The [111] BCC alloys shown in (c-g), along with our predicted hR75-Li19Sn6 not included due to its large size, can be uniquely defined with sequences of Sn blocks, e.g., |3454| for hR48-Li3Sn. These images were generated with the VESTA package [83].
the second shell at 3.20Å (Fig. 5 (b)). Interestingly, previous ab initio studies revealed a lower-energy decoration of the BCC lattice at this composition [47,50]. The simpler hexagonal hP9 phase with the Li 7 Pb 2 prototype does not contain directly interacting Sn atoms and ends up being ∼ 6 meV/atom more stable in (semi)local DFT approximations [47,50,88] and Table S1). Even though the energy difference reduces to ∼ 3 meV/atom at the SCAN level (Table S1), the oS36 phase appears to be only metastable under typical synthesis conditions. We reproduce the previously reported finding [47] that the inclusion of vibrational entropy changes the relative free energy insignificantly (by less than 1 meV/atom up to 800 K in Fig. S2), which is not unexpected given the similarity of the oS36 and hP9 morphologies.
The stable hR33-Li 8 Sn 3 phase [47,50] with the Li 8 Pb 3 prototype has been proposed [47] to explain the nominal compound synthesized in 1996 by Gasior et al. [89]. Just as hP9-Li 7 Sn 2 , it exhibits an ordered pattern along the [111] crystallographic direction of the BCC lattice and can be represented with a hexagonal unit cell (Fig. 5 (d)). In contrast to hP9-Li 7 Sn 2 , it has a peculiar distribution of the Sn atoms: as discussed in Refs. 47 and 50, the occupation of neighboring sites creates Sn-Sn dimers with an unusually short 2.90Å interatomic distance.
The known hP18-Li 13 Sn 5 and hR21-Li 5 Sn 2 phases at nearby compositions bear a strong resemblance to hR33 and contain Sn-Sn dimers of length 2.91Å (on two of the three Sn Wyckoff sites) and 2.92Å (on the only 6c Sn Wyckoff site), respectively. Our PBE calculations agree with the previous results [47] regarding the thermodynamic stability of the former and the metastability of the latter. As summarized in Table S1, hR21-Li 5 Sn 2 is 1.8 meV/atom above the hP18-Li 13 Sn 5 ↔ mP6-LiSn tieline in this approximation but marginally breaks the convex hull by −0.4 meV/atom in the SCAN calculations. The higher sensitivity of the relative stability to the DFT flavor in this case is likely caused by the fairly different morphology of the reference mP6-LiSn phase discussed below.
In order to systematize the description of such BCC structures featuring the same hexagonal base, we introduce a convenient notation that builds on the detailed structural analysis presented in Ref. 47. As can be seen from Fig. 5(c-g), Sn atoms 3, 4, and 5 layers apart along the c axis are shifted by 0, 1/3, and 2/3 unit cell fractions within the a-b base, respectively. Therefore, Sn sequences can be used to uniquely specify these trigonal crystal system structures and easily deduce whether they belong to the hexagonal or rhombohedral lattice system. Indeed, sequences with the total number of Sn layer vertical separations divisible by three (e.g., |45| with 2 Sn atoms in the 9-atom primitive unit cell) result in no net lateral shift and, hence, correspond to hexagonal structures. All other sequences (e.g., |34|) must be tripled to ensure orthogonality of the c axis with the base and, thus, define rhombohedral structures.
At the 3:1 composition, only metastable ambientpressure phases have been proposed so far [47,50] but our evolutionary searches for 12:4 unit cells uncovered a new thermodynamically stable member of the BCC family, hR48-Li 3 Sn with the |3435| sequence. The rhombohedral structure with a 16-atom unit cell appeared in both NN evolutionary runs after 849, and 1,491 local op- timizations, although it placed 17 meV/atom above a cF16 BCC phase and ranked 4th in the NN pool. In the considered DFT approximations, it was found to be the lowest-energy 3:1 configuration and stable with respect to the decomposition into the previously proposed hP9-Li 7 Sn 2 and hR33-Li 8 Sn 3 BCC-based ground states by about −2 meV/atom at T = 0 K (Table S1). Following phonon calculations indicated that hR48-Li 3 Sn should remain stable at elevated temperatures. The finding motivated us to examine BCC structures more closely for other possible ground states. The cluster expansion method has been widely used to explore combinatorially large configuration spaces on given lattices [90,91] but an exhaustive screening of all possible decorations of large unit cells would require a separate dedicated study. Fortunately, the observed favorability of the [111] BCC motifs allowed us to focus on the most promising subspace. Instead of relying on traditional atomic clusters, we expanded the total energy as a linear combination of various blocks comprised of 3, 4, or 5 layered units described above: where c i is the energy of block i and N i is the number of such blocks in a structure.  (Fig. 6). Based on this information, we constructed a number of new members with up to 6 units and plotted their DFT relative formation energies (Fig. 7). This analysis led to identification of yet another thermodynamically stable BCC phase, hR75-Li 19 Sn 6 . Comprised of 6 layered units, the |345454| structure is far more difficult to find with ab initio searches because it has an unusual composition and a 25-atom primitive rhombohedral unit cell with a large c/a ≈ 15 ratio. Our three NN evolutionary runs with 10 4 local relaxations failed to locate the presumed 19:6 global minimum and converged to a monoclinic mS50 structure 2.1 meV/atom (7.3 meV/atom in the PBE) above hR75. According to our DFT tests (Table S1), hR75-Li 19 Sn 6 lies essentially on the hP9-Li 7 Sn 2 ↔ hR48-Li 3 Sn tie-line and below the hP9-Li 7 Sn 2 ↔ hR33-Li 8 Sn 3 tie-line by 1-2 meV/atom. Alloys with higher Sn content (x ≥ 0.3) adopt stable configurations with distinctive extended Sn fragments or frameworks [50]. The observed mP20-Li 7 Sn 3 phase [92] is a peculiar decoration of the BCC lattice resulting in Sn zig-zag trimers with relatively short bonds of 2.99 and 3.00Å. It has been found slightly metastable by 1-2 meV/atom with respect to hP18-Li 13 Sn 5 and mP6-LiSn in the standard DFT calculations (Table S1) [47,50]. Our phonon calculations suggest a stabilizing effect of the vibrational entropy but the flatness of the relative free energy curve makes it difficult to give a reliable estimate of the phase transition temperature (Fig. 3). Moreover, our SCAN calculations show that it may be a true T = 0 K ground state.
The previously proposed oS12-Li 2 Sn phase [50] was found to be the most stable 2:1 configuration in our calculations as well. The structure contains planar zig-zag Sn chains separated by two Li layers. In contrast to other Li-Sn phases discussed so far, the Li and Sn sites have 12-atom coordinations and are distributed on an HCP lattice. The PBE calculations with an ultrasoft pseudopotential placed it just 1 meV/atom above the hP18-Li 13 Sn 5 ↔ mP6-LiSn tie-line [50]. Our PBE and SCAN results (Table S1) indicate that oS12-Li 2 Sn may actually be thermodynamically stable at zero and elevated temperatures (Fig. 3).
The experimental mP6-LiSn phase [93] is another example of a stable BCC alloy. A clustered population of the lattice sites creates fully connected zig-zag Sn nets shown in Fig. 5(j). Our evolutionary searches for various unit cell sizes of the 1:1 composition located several DFT local minima within a 20 meV/atom energy window (Figs. 3 and S1). The most dramatic entropy-driven stabilization relative to mP6, from 14.9 meV/atom at 0 K down to 6.4 meV/atom at 600 K, occurred for a tI24 phase observed in Ref. 94, but none of the considered polymorphs became stable in our calculations (Figs. S1 and S2) up to the reported melting temperature of 758 K at this composition [95].
Our screening of competitive high-T compounds produced a viable FCC-based mS40-LiSn 4 phase (Fig. 5(l)) that breaks the tP14-Li 2 Sn 5 ↔ α-Sn tie-line at 510 K. This temperature estimate is particularly sensitive to the description of pure Sn. The difficulty of reproducing α → β transition temperatures and pressures for group-XIV Si and Sn has been pointed out in several studies [70][71][72][97][98][99][100]. The overestimated temperature value for the α-Sn to β-Sn transformation obtained in our standard PBE calculations is 540 K. A lower free energy of the reference β-Sn could shift the mS40-LiSn 4 stabilization to significantly higher temperatures.
To the best of our knowledge, there have been no attempts yet to synthesize Li-Sn alloys at high pressures, which makes Sen and Johari's ab initio study [47] an important baseline for stability of binary compounds up to 20 GPa. Our evolutionary searches performed at 20 GPa reproduced all the previously reported ground states for compositions above x = 2/9. The 7:2, 8:3, 13:5, and 5:2 compounds retain their ambient-pressure BCC prototypes discussed above, while 3:1 and 1:1 compounds adopt well-known simpler cF16 (D0 3 ) and cP2 (B2) prototypes, respectively. Compared to the ambient-pressure results, there are very few competing phases within the 20 meV/atom window above x = 1/3. The closest to stability phase is tI12-LiSn 2 , which is 1 meV/atom above cP2-LiSn ↔ BCT-Sn as seen in Fig. 4.
For compositions below x = 2/9, our searches produced several new stable Li-Sn phases that reshape the Li-rich part of the previously proposed convex hull at 20 GPa. Given the relatively low symmetry of these phases and the apparent lack of well-defined underlying motifs, their structural features are better illustrated with the RDF plots shown in Fig. S3. It can be seen that Sn atoms in these alloys are separated by at least 4.0Å, and the Sn local environments consist of 14-17 Li atoms 2.48-3.06Å away, with more Sn-rich phases generally having shorter Sn-Li interatomic distances. Despite the morphological commonalities, the slight variations lead to noticeable differences in enthalpies. At the 7:1 stoichiometry, we obtained several competing configurations with low symmetry. The most stable aP32 structure is −1.9 meV/atom below aP16 proposed previously [47]. Both aP16 and aP32 have Sn sites with well-defined shells of 16 Li atoms within the 2.6-3.0Å range. By linking Sn atoms with 4.5-5.0-Å connections, one can discern different extended Sn frameworks with channels of variable size along a, larger in aP32, filled with Li atoms.
According to our global optimization results, the most Li-rich ground state at 20 GPa actually occurs at the 11:2 composition and has a low-symmetry aP26 structure. It makes all 7:1 phases metastable by at least 3.7 meV/atom. While both the NN and DFT methods favor the aP26 polymorph, the latter indicates that another 11:2 phase, mP26, is only 0.4 meV/atom above the putative ground state at T = 0 K and becomes stable at 460 K. Notwithstanding the near degeneracy in the full temperature range, the two structures feature different Li coordinations of the two Sn sites: 16 and 16 in aP16 versus 15 and 17 in mP26.
The most significant improvement in stability with respect to the previously proposed ground states is observed at the 5:1 stoichiometry. The best configurations identified in previous [47] and our studies have the same Pearson symbol (mS48) and space group (C2/m) but the latter is more stable by 10.7 meV/atom; in fact, the former is only metastable by 8.9 meV/atom relative to the new set of ground states at neighboring compositions (Table S2). Morphologically, the previous phase features only 4i Li/Sn Wyckoff positions, Sn sites coordinated with 15 and 17 Li atoms 2.61-3.06Å away, and a Sn framework with Li-filled channels along b, whereas ours contains four 8j Li Wyckoff positions, Sn sites coordinated with 14 Li atoms only 2.55-2.75Å away, and a Sn framework without apparent extended pores.
A mS44-Li 9 Sn 2 phase found in our study represents another stable compound at a new stoichiometry, being at least −8 meV/atom lower in enthalpy with respect to the best 5:1 and 15:4 phases. Compared to the compressed structures discussed so far, mS44 has Sn-Li interatomic distances dispersed over a wider 2.50-3.36Å range. The inclusion of this phase in the convex hull destabilizes the previously proposed tI50-Li 4 Sn [47] by about 1 meV/atom (see Table S2).
Finally, we sampled the unusual 15:4 composition because phases with the cI76-Cu 15 Si 4 prototype (I43d) have been observed in related Li-Si and Li-Ge binaries [101,102] and because another oS152-Li 15 Si 4 prototype (F dd2) has been shown to become stable in the Li-Si system above 7.9 GPa [103]. Our DFT calculations at 20 GPa indicate that the cubic structure is unstable by over 60 meV/atom but the orthorhombic one is actually −1.3 meV/atom below the tie-line connecting previously proposed 5:1 and 7:2 ground states [47] and could have been considered stable. We find that oS152-Li 15 Sn 4 breaks the new 9:2 and 7:2 tie-line as well, albeit by only −0.6 meV/atom. However, our evolutionary searches identified a new hR57 structure −1.8 meV/atom lower in enthalpy compared to oS152. This presumed T = 0 K ground state with a large rhombohedral unit cell belongs to the [111] BCC family and can be represented simply as |4555|. Our phonon calculations demonstrate that oS152 does stabilize over hR57 above 400 K.
The impressive number of viable phases with large lowsymmetry unit cells appearing at ambient and high pressures testifies to the exceptional complexity of the Li-Sn phase diagram. While we find it encouraging that our extensive screening accelerated with machine learning uncovered several new phases with lower Gibbs free energies, it cannot be considered exhaustive. Future explorations of the Li-Sn system with higher-accuracy MLPs and DFT approximations may lead to further revision of the phase diagram.

C. Analysis of electronic and other properties of ambient-pressure Li-Sn compounds
Electronic, mechanical, and electrochemical properties of Li-Sn alloys have been investigated in several ab initio studies [47][48][49][50]104]. Analyses of the projected DOS and Bader decomposition data have indicated the expected metallicity of all Li-Sn alloys and transfer of charge from Li to Sn. Interestingly, the charge redistribution in alloys with high Li concentrations exceeding 4:1 has been shown to leave Li atoms in either ionized (+0.8) or nearly neutral (−0.3 -0.0) states [47]. The manifold of the electronic states near the Fermi level has primarily the Sn-p character but displays increasing Li-s/p contributions in Li-rich compounds [47,49]. It has been pointed out that addition of Li makes the dominant bonding in Li-Sn alloys evolve from covalent Sn-Sn/Li-Sn with a characteristic pseudo-gap to metallic Li-Li [49].
Examination of DOS profiles near the Fermi level can indeed provide insights into materials' various properties ranging from transport, superconductivity, and mag- netism to stability. The presence of a pseudo-gap and the population of bonding states has been found to correlate with structure's favorability in numerous cases [5,76,105]. For instance, our studies of related chemical systems have shown that the unusual stoichiometry of LiB x (x ∼ 0.9) and the stability of NaSn 2 result from the placement of the Fermi level near the bottom of a well-defined DOS valley [5,76,106]. Figure 8 shows that most of the known and proposed Li-Sn alloys, especially those with x ≥ 0.5, have a significant DOS at the Fermi level and may not conform to this simple model. In fact, it is not easy to disentangle the hybridization and charge transfer factors that determine the shape and the population of the electronic states with the Sn-p and Li-s/p character. The Sn-s states, on the other hand, exhibit little hybridization and provide information about the behavior of Sn-centered orbitals. With the large s-p splitting disfavoring the formation of the traditional sp 2 /sp 3 hybrids in pure Sn open structures, the low-lying s states split first and only then does the antibonding s set mixes partially with the p states ( Fig.  8(o)) [107,108].
In the β-Sn, mS40-LiSn 4 , tP14-Li 2 Sn 5 , and mP6-LiSn phases, the Sn-s states disperse into a single band that partially overlaps in energy with the bottom of the Snp band. As the number of Sn nearest atoms within the 3.08-3.40Å range changes from 6 to 8-12, 6, and 4 in this series, the Sn-s bandwidth gradually decreases. In oS12-Li 2 Sn with zig-zag chains, the dispersion reduces down to about 4 eV which creates a 1-eV separation between the Sn-s set and the Sn-p and Li-s/p manifold. All shown phases with higher Li content do not have extended Sn frameworks and feature fairly localized Sn-s states. The number and the width of the DOS peaks are determined by the distribution of isolated Sn trimers (in mP20-Li 7 Sn 3 ), Sn dimers (for BCC phases with 0.24 < x < 0.277), and single Sn atoms (in hP9-Li 7 Sn 2 and cF420-Li 17 Sn 4 ). For example, the presence of the 2.9-Å Sn-Sn dimers within the type-3 units in the [111]-BCC phases (Fig 5(e-h)) produces two 0.7-eV peaks centered about 2 eV apart, while the orbital overlaps between Sn atoms 4.7Å apart within the [111]-BCC base or the type-4 blocks give rise to the central 1-eV peaks (Fig 8(c,e-h)). The lower edge of the Sn-p band moves in concert with the bandwidths of the Sn-s states, shifting from −4 eV in oS12-Li 2 Sn to −3 eV in cF420-Li 17 Sn 4 . Interestingly, the lowest states in the manifold in oS12-Li 2 Sn and hP9-Li 7 Sn 2 have the Li-s/p character due to a larger dispersion of these nearly free electron states that can be seen more clearly on band structure plots for select phases (Figs. S4-S6).
The change in volume upon Li insertion is one of the critical factors for durability of Li-ion battery anodes. It has been shown that the calculated atomic volume in Li-Sn alloys notably deviates from the linear dependence on the Li content [47,48]. Our results in Fig. 9 demonstrate that the atomic volume correlates with the formation energies across the Li 1−x Sn x composition range. The lowest values occur in compounds with the maximum transfer of charge from Li to Sn near the 4:1 stoichiometry.
A similar trend in the full composition range can be observed for vibrational properties using the phonon DOS results in Fig. 8 Li 2 Sn 5 based on a 3D Sn framework, and eventually settles at 12 meV in BCC-based Li-rich alloys. The highest frequencies of the Li modes actually increase with the first addition of Sn, reaching a maximum of 51 meV in cF420-Li 17 Sn 4 , and then drop precipitously in alloys with x ≥ 0.5. Plotted versus the atomic volume in Fig. 9(b), the highest phonon frequencies in all Li-Sn phases exhibit a clearer trend: the BCC-Li and α-Sn outliers fall closer to the linear fit, which leaves tP14-Li 2 Sn 5 with hard Li modes as the only significant deviation.
Finally, we comment on how the predicted ambientpressure phases could be detected. Electrochemical cycling measurements have been used to identify Li-Sn alloys in several studies [109][110][111]. We update the previously reported average voltage plot [50] by including the hR75-Li 19 Sn 6 and hR48-Li 3 Sn phases predicted in our study. Since these phases deepen the convex hull by only 1-2 meV/atom, the additional features near y = 3 in Li y Sn (Fig. 9) are very fine and would be difficult to observe. However, their powder x-ray diffraction (XRD) patterns differ considerably from those in previously observed or proposed BCC phases and would provide strong evidence of their formation (Fig. S7). The combination of electrochemical and XRD measurements have been recently used to detect a new NaSn 2 phase [112] predicted in our earlier ab initio study [76].

IV. SUMMARY
Our re-examination of the Li-Sn binary illustrates that well-known systems may still host new synthesizable compounds and that computationally demanding ab initio prediction of thermodynamically stable materials can be effectively boosted with MLPs. Given the surprising scarcity of successful MLP-assisted predictions of stable compounds so far, our study was focused on devising guidelines for constructing reliable models, identifying viable structures, and establishing stability trends that increase the chance of finding ab initio ground states. To this extent, the combination of the evolutionary sampling and stratified training allowed us to build a practical Li-Sn NN potential with a relatively modest 10.2 meV/atom accuracy but a reliable performance in unconstrained searches. The model construction involved less than 2% of the total computational cost required for the systematic exploration of the binary system. Our screening of over 1.1 million (2.3 thousand) crystalline phases at zero (elevated) temperature with the NN potential along with the following DFT analysis of select candidates cost only an estimated 0.1-1% of the traditional ab initio searches. As a result, we discovered overlooked thermodynamically stable hR48-Li 3 Sn and mS40-LiSn 4 at ambient pressure and several possible new ground states at 20 GPa. In contrast to the previously identified stable hP9-Li 7 Sn 2 and hR33-Li 8 Sn 3 phases [47,50] that have known Li 7 Ge 2 and Li 8 Pb 3 prototypes in related binaries and can be naturally identified by chemical substitution, our large hR48 and mS40 structures appear to have no analogues in other chemical systems. Moreover, after introducing a phenomenological model for the [111] BCC family of Li-rich alloys, we constructed another stable hR75-Li 19 Sn 6 phase with an even larger unit cell that could not be found with our standard evolutionary searches. These findings highlight the benefit of combining different prediction strategies that utilize available experimental knowledge, explore configuration spaces with global optimization algorithms, and involve rational design based on stability analysis. Since the thermodynamic stability of the new compounds predicted with this method is ultimately determined by the accuracy of the employed DFT approximation, we performed additional tests with various DFT flavors. Our results indicate that the relatively small formation enthalpy differences of a few meV/atom are fairly insensitive to the systematic and numerical errors due to the similarity of the competing configurations. For the same reason, the vibrational entropy contribution to the Gibbs free energy was found to be an insignificant factor for differentiating Li-rich phases with related BCC motifs. Consequently, the vast majority of detected metastable phases in the 20 meV/atom enthalpy window, a common criterion for viable candidates, proved to be irrelevant in the Li-Sn system. The proposed thermodynamically stable phases are expected to be synthesizable and detectable with standard powder XRD measurements. The formation of hR48-Li 3 Sn and hR75-Li 19 Sn 6 would lead to further refinement of the electrochemical profile at the Li-rich end of the complex phase diagram.      . Powder x-ray diffraction patterns simulated for λCu K−α = 1.54059Å for some of the known and predicted ambientpressure Li-Sn phases.