Discovery of oxide Li-conducting electrolytes with corner-sharing frameworks via topology-constrained crystal structure prediction

Abstract

Oxide Li-conducting solid-state electrolytes (SSEs) offer excellent chemical and thermal stability but typically exhibit lower ionic conductivity than sulfides and chlorides. This motivates the search for new oxide materials with enhanced conductivity. Crystal structure prediction is a powerful approach for identifying such candidates. However, the structural complexity of oxide SSEs, often involving unit cells with more than 100 atoms, presents significant challenges for conventional methods. In this study, we introduce TOPIC, a structure prediction algorithm that reduces configurational complexity by enforcing corner-sharing (CS) bond topology constraints. We demonstrate that TOPIC successfully reproduces the ground-state and metastable structures of known oxide SSEs, including LiTa₂PO₈ and Li₇La₃Zr₂O₁₂, which contain up to about 200 atoms per unit cell. By combining this approach with a pretrained machine-learning interatomic potential, we systematically screen quaternary oxide compositions and identify 92 potential candidates with CS frameworks. In particular, Li₄Hf₂Si₃O₁₂, which corresponds to the ground state at its composition, exhibits an ionic conductivity of 14 mS cm⁻¹ (7.56–26.22 mS cm⁻¹, considering uncertainty), a hull energy of 21 meV atom⁻¹, and a band gap of 6.5 eV. Through our investigation, we identify the Li ratio as one of the key factors determining the stability of CS structures. Overall, our approach provides a practical and scalable pathway for discovering high-performance oxide solid electrolytes in previously unexplored chemical spaces.

---

1 Introduction

Despite extensive research into lithium-ion battery technologies, persistent safety issues and limited capacity continue to impede their widespread adoption in high-energy applications, such as electric vehicles and large-scale energy storage.^1–5 Solid-state electrolytes (SSEs) offer an attractive solution by supporting stable operation at higher voltages.^6,7 Among SSE candidates, oxide-based electrolytes are particularly promising due to their superior chemical and thermal stability.⁸ Several oxide SSEs, such as garnet-type Li₇La₃Zr₂O₁₂ (LLZO),⁹ perovskite-type Li_3xLa_(2/3−x)TiO₃ (LLTO),¹⁰ NASICON-type Li_1.3Al_0.3Ti_1.7(PO₄)₃ (LATP),¹¹ and LiTa₂PO₈ (ref. 12), have achieved ionic conductivities above 1 mS cm⁻¹. However, their ionic conductivities are still generally lower than those of sulfide materials,¹ highlighting the need to discover novel oxide SSEs.

Computational screening can accelerate discovery of new SSE materials at a much higher speed compared to experiments. In fact, there have been efforts to discover SSEs by searching existing databases, such as the Inorganic Crystal Structure Database (ICSD),¹³ which have successfully identified promising candidates.^14–16 To further expand the search space, several studies have employed structures available in existing databases as templates, systematically substituting their constituent elements to generate new candidate materials.^17–20 However, these data-driven schemes inherently lack the capability to uncover new structural prototypes, posing a critical challenge for exploring uncharted chemistries.²¹

Discovery of materials with novel structural prototypes requires crystal structure prediction, which is a computational approach used to identify ground-state or metastable crystal structures for a given chemical composition.^22–24 These methods have successfully led to the discovery of diverse material applications, such as the superhard materials,²⁵ high-temperature superconductors,²⁶ and two-dimensional electrides.²⁷ Crystal structure prediction is also used to discover novel SSE materials.^28–31 However, the discovery of oxide SSEs via crystal structure prediction remains significantly more challenging than for other types, such as sulfides or chlorides. This is because oxide SSEs often possess highly complex structures, typically containing over 100 atoms per unit cell and involving quaternary or higher-order compositions.¹ Conventional structure prediction methods, including genetic algorithms, are generally restricted to simpler systems because of the high computational cost for evaluating energies of intermediate and candidate structures using density functional theory (DFT) calculations. To address this limitation, several studies have employed machine-learning interatomic potentials (MLIPs)^32–34—surrogate models trained on DFT data to predict energies and forces—in crystal structure prediction, as they offer significantly higher speed and comparable accuracy to DFT calculations.^35–38 In fact, MLIP-accelerated crystal structure prediction has recently been explored for the discovery of novel SSE materials.³⁹ However, even when accelerated by MLIPs, reported crystal structure prediction studies have thus far been limited to ternary systems with fewer than 50 atoms per unit cell, or to multinary systems with high-symmetry structures,⁴⁰ which remain insufficient for the discovery of oxide SSE materials. Therefore, more efficient prediction methods are needed to address the complexity of oxide SSE structures.

Using structural constraints in crystal structure prediction algorithms can further accelerate the process, especially when prior knowledge of structural characteristics (such as the shapes of known polyanions) is available^30,41 or when the search space can be adaptively reduced based on symmetry principles.⁴² In the case of oxide SSEs, oxygen-corner sharing (CS) frameworks have already been recognized as a promising structural motif for fast Li-ion conduction,^12,43 whereas sulfides and halides typically include isolated units such as PS₄. Based on this knowledge, experimental efforts targeting compositions likely to yield CS frameworks led to the discovery of LiTa₂PO₈, which was indeed confirmed to exhibit a CS framework.¹² Furthermore, computational screening of the Materials Project database⁴⁴ focusing on CS-framework structures successfully identified oxide materials with high ionic conductivity, including LiGa(SeO₃)₂, which was also validated through experiments.⁴³ In that study, the superior ionic conductivity of CS frameworks was rationalized through theoretical analysis, providing insight into their structure–property relationships. However, these conventional theoretical prediction methods are still limited in their ability to discover oxide materials because of the large structural size and low symmetry of these systems.

Here, motivated by the previously established importance of CS frameworks in oxide lithium-ion conductors, we develop a crystal structure prediction algorithm named TOPIC (TOpology-constrained Prediction of Inorganic Crystals). TOPIC generates candidate structures under CS bond topology constraints, significantly reducing the configurational search space compared to conventional methods. We first validate the TOPIC algorithm by successfully reproducing the crystal structures of known oxide systems containing up to 192 atoms per unit cell, including LLZO. Next, using a pretrained MLIP (SevenNet-0),⁴⁵ we systematically screen for novel quaternary oxide SSEs exhibiting CS frameworks generated by TOPIC. Our initial screening focuses on chemical systems composed of commonly used elements, such as Ti, Zr, and P, where we find that the lithium content serves as a key descriptor for the stability of CS frameworks. Based on this descriptor, we enumerate possible quaternary compositions likely to exhibit CS frameworks composed of octahedral and tetrahedral units and apply TOPIC to explore such chemical space. Overall, we identify 92 new candidate materials with potentially high ionic conductivity. Among them, 15 oxide materials are further validated using DFT calculations to assess their thermodynamic stability, electronic band gaps, and ionic conductivity. Finally, we report several novel framework types with high predicted ionic conductivity and derive design principles by analyzing the predicted SSE structures and the Li-ion conduction pathways.

2 Results

2.1 Definition of corner-sharing bond topology

We define a bond topology as the connectivity of a three-dimensional network formed by polyhedra centered on cations (except Li) with surrounding oxygen atoms. A fully connected CS framework is one in which every oxygen bridges exactly two polyhedra, thereby establishing continuous connectivity (left panel of Fig. 1a). This definition explicitly excludes configurations where an oxygen anion binds to only one polyhedron, forming an isolated vertex (middle panel of Fig. 1a), or configurations where two polyhedra share an edge or face (right panel of Fig. 1a). Our definition thus slightly differs from the broader criterion proposed by Jun et al.,⁴³ which allows isolated vertices. Hereafter, we refer to frameworks with fully connected vertex-sharing topologies as CS frameworks, and we classify all other cases as non-CS frameworks.

Fig. 1 TOPIC algorithm. (a) Schematic illustrations of a fully connected corner-sharing framework (left), a framework with an isolated vertex (middle), and a framework with edge-sharing polyhedra (right). The dashed red circle in the middle panel indicates the isolated vertex in the framework, and dashed green circles in the right panel indicate nodes on edge-shared parts. Oc and T represent the cations forming octahedral and tetrahedral polyhedra, respectively. (b) Overview of the TOPIC algorithm; the LiTiPO₅ case is shown as an example. (c) Detailed schematic process of cation–oxygen framework generation. Spring lines in the second panel indicate the virtual bonds between cation–oxygen pairs, to which the Lennard-Jones potential is applied.

2.2 TOPIC algorithm

We develop TOPIC, a random structure generation algorithm constrained by CS bond topology. Structural optimizations and energy evaluations are carried out using MLIPs during the structure search, and the final candidates are recalculated with DFT to obtain accurate energies. To construct the training set for MLIPs, we perform melt–quench–annealing molecular dynamics (MD) simulations at a target SSE composition to generate disordered structures using the SPINNER code,⁴⁰ following a procedure similar to our previous studies.^24,36 From this training set, we train three types of MLIPs for each composition: (1) a 3 Å cutoff model excluding Li atoms, (2) a 6 Å cutoff model excluding Li atoms, and (3) a 6 Å cutoff model including all atoms. The first two models are trained by removing Li atoms from the configurations while retaining the original energy and force labels for the others, which include the contribution from Li. In other words, the energetic contribution of Li is implicitly averaged out during training. These three MLIPs are then used in the structural search. Interatomic distances (such as cation–oxygen bond lengths) are extracted from the MD trajectories for generating training set and used as geometric constraints during structure generation (see Methods for details).

TOPIC generates random CS frameworks without Li and subsequently determines the optimized Li positions. Specifically, the algorithm sequentially performs the following steps: (1) generating cation sublattices, (2) generating the cation–O framework by placing oxygen atoms under bond topology constraints followed by relaxations, and (3) placing Li atoms (see Fig. 1b). At the first step, the cation sublattices are generated using random structure generation with space group constraints using the RandSpg code,⁴⁶ as similar methods are commonly adopted in other structure prediction algorithms.^22,47,48 At the final step, Li atoms are placed via Monte Carlo (MC) simulations coupled with MLIP-based optimizations. In this approach, virtual bonds are first formed between selected cation pairs on the given cation sublattice, and oxygen atoms are placed at the midpoints of these bonds (see Fig. 1c). The resulting network is required to satisfy the CS bond topology and the predefined coordination number for each cation species, as determined by Pauling's rules⁴⁹ and Shannon's ionic radii⁵⁰ (Table S1). (We also confirm that O atoms are located near the midpoint between cation pairs in experimentally known materials; see Fig. S1). The resulting structures are optimized by applying a Lennard-Jones (LJ) potential to the virtual bonds to enforce reasonable cation–oxygen bond lengths. In addition, a harmonic repulsive potential is used to prevent unphysically short distances and unintended increases in coordination numbers during optimizations. Finally, the structures are further relaxed using MLIP optimizations. After each relaxation step (both LJ and MLIP), we verify whether the structure still satisfies the CS bond topology and the predefined coordination numbers. Structures that violate these criteria are discarded. Further details of each algorithmic step are provided in the Methods section.

In the CS-framework-generation step, structural optimizations are performed using MLIPs trained on datasets that exclude the effects of Li atoms, as described above. With a cutoff radius of 6 Å, these MLIPs yield average validation errors of 0.014 eV atom⁻¹ for energy and 0.87 eV Å⁻¹ for force. More accurate calculations are subsequently conducted with MLIPs trained including Li, during the Li-insertion stage. This two-step optimization procedure (from Li-free to Li-occupied systems) effectively identifies structures that satisfy CS-topology constraints while reducing computational cost, as demonstrated in the following subsection. The effectiveness of the initial framework-generation step arises from the intrinsically weak Li–O interactions in multicomponent Li oxides, which result in only minor modifications of the Li-free structural framework upon Li insertion. Weak Li–O interactions have been reported in NASICON-type Li compounds, which exhibit low-frequency phonon modes and notably weak Li–O bonding, as confirmed by crystal orbital Hamilton population (COHP) analyses.^51,52 Our COHP analysis similarly confirms weak Li–O interactions in compounds such as LiTi₂(PO₄)₃, LiGa(SeO₃)₂, and LiTa₂PO₈, as well as in the melt-quenched amorphous phases of the same compositions (Fig. S2).

2.3 Validation of TOPIC

To assess the reliability of TOPIC, we check whether it reproduces the structures of known materials. The test set includes three experimentally reported compositons (LiTi₂(PO₄)₃, LiAlSiO₄, and LiTa₂PO₈); four computationally proposed candidates (LiTaSiO₅, LiTiPO₅, Li₂Mg₂(SO₄)₃, and LiGa(SeO₃)₂);^43,53 and one garnet-type Na compound with corner-sharing bond topology (Na₃Ga₃Te₂O₁₂). For the two compositions, LiTiPO₅ and LiAlSiO₄, which exhibit multiple experimentally observed phases, we examine whether TOPIC accurately reproduces all of these phases. Specifically, LiTiPO₅ exhibits α (P1̄) and β (Pnma) phases,⁵⁴ while LiAlSiO₄ has three known phases: α (R3), β (P6₄22), and γ (Pc).⁵⁵ Therefore, total eleven structures are accounted for in the test. We note that LiTa₂PO₈ exhibits disordered Li sites, whereas the other compounds have fixed Li sites. A summary of test systems and search outcomes is provided in Table 1.

Table 1 Information of validation test systems. The space group, the number of atoms in the unit cell (N_atom), and the number of discovered reference structures in 300 000 trials (N_found) are provided

Formula	Space group	Natom	Nfound
LiTi2(PO4)3	R3̄c (167)	108	7620
LiGa(SeO3)2	I4̄ 2d (122)	80	24
Li2Mg2(SO4)3	Pbcn (60)	76	2
LiTa2PO8	C2/c (15)	96	2
Na3Ga3Te2O12	Ia3̄d (230)	160	10 575
α-LiTiPO5	P1̄ (2)	32	1
β-LiTiPO5	Pnma (62)	32	155
α-LiAlSiO4	R3 (146)	126	148
β-LiAlSiO4	P6422 (181)	84	12
γ-LiAlSiO4	Pc (7)	28	1
LiTaSiO5	P21/c (14)	32	76
Li7La3Zr2O12	Ia3̄d (230)	192	10 512

---

For each system, we perform 300 000 trials and identify the lowest-MLIP-energy structure, followed by DFT evaluations of candidates within a 30 meV atom⁻¹ window of this lowest-energy structure. Note that the number 300 000 refers to the structures initially generated to satisfy the CS bond topology, prior to any optimization with the LJ potential or MLIPs. After optimization, only a small fraction of these 300 000 structures remain (approximately 7%), preserving the CS framework (see Methods for details).

Fig. 2a presents the lowest-energy structures (within DFT calculations) obtained for each system, and comparisons with ICSD reference structures (including metastable polymorphs) are shown in Fig. S3. Among ten target structures with ordered Li sites, TOPIC successfully predicts six with the correct Li positions: LiTi₂(PO₄)₃, LiGa(SeO₃)₂, β-LiTiPO₅, α-LiAlSiO₄, β-LiAlSiO₄, and Na₃Ga₃Te₂O₁₂. For Li₂Mg₂(SO₄)₃, the predicted framework reproduces the known structure with only minor deviations in Li sites, yet all predicted Li positions remain within established conduction pathways (Fig. S3g). For LiTaSiO₅, α-LiTiPO₅, and γ-LiAlSiO₄, TOPIC predicts slightly distorted structures, while the polyhedral connectivities remain consistent with the reference structures. In the case of LiTa₂PO₈, the reference structure contains partially occupied Li sites; here, we confirm that the predicted framework is identical to the reference and that the predicted Li positions are on the conduction channels (Fig. S3h).

Fig. 2 Validation of the TOPIC algorithm. (a) The lowest-energy structures for each validation test system generated by TOPIC. (b) The lowest-energy and a new metastable structure in the LiTa₂PO₈ system obtained by TOPIC. (c) Diffusion coefficients at 400, 700, and 1000 K for all polymorphs of the LiTa₂PO₈ system obtained by TOPIC. Blue lines represent the known structure, the red line corresponds to the structure of the potential candidate, and gray lines indicate other polymorphs. Solid lines denote results calculated by MD simulations using SevenNet-0, while dashed lines denote results obtained by DFT simulations. (d) Revised structure generation process for cubic Li₇La₃Zr₂O₁₂. Left: the framework consisting of ZrO₆ octahedra and LiO₄ tetrahedra generated by TOPIC. Middle: the framework after La atoms are introduced via Voronoi tessellation and Monte Carlo simulation. Right: the final Li₇La₃Zr₂O₁₂ structure after the Li insertion.

It is worth noting that TOPIC enables the more efficient identification of complex crystal structures that are challenging to discover using conventional evolutionary algorithms, which are widely adopted for crystal structure prediction. To assess the practical advantage of TOPIC, we benchmark it against SPINNER,^24,36,40 an MLIP-accelerated evolutionary algorithm without topology constraints, using the same number of generated structures and comparable CSP settings (SI Note 1 and Table S2 for details). For all benchmark systems examined, TOPIC discovers the reference structures more frequently and with shorter wall time than SPINNER. We also compare TOPIC with MAGUS, a state-of-the-art structure prediction method based on symmetry principles.⁴² We select Mg₃Al₂Si₃O₁₂ as a benchmark system, as it was well examined in ref. 42 and adopts a corner-sharing framework when Mg atoms are excluded. TOPIC found the ground-state structure approximately once every 28 generated structures, whereas MAGUS found it approximately once every 125 generated structures. These results demonstrate the clear advantage of TOPIC's topology-constrained search strategy and highlight the importance of incorporating topology constraints for the objectives of this study.

TOPIC also predicts low-energy structures different from those previously reported. For example, in LiTaSiO₅, it identifies a structure that is 6 meV atom⁻¹ lower in energy than the reported phase, distinguished by a different stacking sequence (Fig. S4a). In addition, we discover a previously unreported metastable polymorph of LiTa₂PO₈, which lies 17 meV atom⁻¹ higher in energy than the reported phase. This polymorph differs from the known phase in its stacking sequence but retains similar polyhedral motifs (Fig. 2b and S4b). Note that this new polymorph of LiTa₂PO₈ exhibits nearly twice the conductivity of the previously known phase (Fig. 2c, dashed line).

Finally, we discuss LLZO, which is the only oxide SSE known not to exhibit a CS framework. However, if one considers a framework built from ZrO₆ octahedra and LiO₄ tetrahedra instead of focusing on Zr and La, the cubic phase of LLZO can be interpreted as a corner-sharing framework (i.e., Li₄La₃(ZrO_6/2)₂(LiO_4/2)₃).¹² We slightly modify the TOPIC algorithm (see Methods for details) to examine whether it can reproduce the structure of LLZO (see Fig. 2d). Specifically, we first predict the ZrO₆–LiO₄ framework, and then place La atoms using MC simulations. Finally, we remove the pre-existing Li atoms and re-predict all Li positions with MC simulations. From this modified process, we successfully identify the LLZO structure, thereby confirming the effectiveness of the TOPIC algorithm. In the screening process outlined in the following subsections, La is excluded from the search space because the method requires modification. Future research, however, could systematically extend the TOPIC algorithm to this material family, as demonstrated above.

2.4 Structure exploration in representative Li-SSE element sets

We first explore candidate materials composed of cation polyhedra with octahedral and tetrahedral coordination, with general stoichiometry Li_x(OcO_6/2)_m(TO_4/2)_n. Here, the lithium content (x) is determined by charge neutrality, and Oc and T denote cations occupying octahedral and tetrahedral sites, respectively. Based on their frequent occurrence in oxide SSE frameworks,^8,56 we select Ta, Ti, Zr, and Al for Oc, and P, Si, and Al for T. Considering Oc : T ratios of 2 : 1, 3 : 2, 1 : 1, 2 : 3, and 1 : 2—the most common cation ratios observed in Li-containing quaternary oxides (see Fig. S5)—we generate 50 candidate compositions in total. From this set, we direct our attention to new compositions that have not yet been explored experimentally and are absent from ICSD.¹³ After removing compositions that cannot satisfy charge neutrality for x > 0, 44 unique compositions remain as our target space. For comparison, we apply SPINNER to generate structures beyond the CS topology. SPINNER generates 60 000 candidate structures at Z = 4 for Oc : T ratios of 2 : 1, 1 : 1, and 1 : 2, and Z = 2 for ratios of 2 : 3 and 3 : 2, resulting in unit cells with approximately 50 atoms to maintain computational feasibility. Following structural generation in both methods, DFT energies are computed for candidates lying within 50 meV atom⁻¹ of the lowest-MLIP-energy structure.

By comparing the lowest-energy structures obtained from each method, we calculate the energies above the convex hull E^DFT_hull for both CS and non-CS frameworks, hereafter denoted as E_hull^DFT (CS) and E_hull^DFT (non-CS), respectively. Fig. 3a shows E_hull^DFT (CS) and E_hull^DFT (non-CS) across all compositions studied in the present work. Note that experimentally synthesized SSEs often exhibit positive E_hull^DFT values (e.g., LiTa₂PO₈: 26 meV atom⁻¹).^57,58 In comparison, the SSEs identified in this work also exhibit positive but comparable or smaller E_hull^DFT values, suggesting their potential synthesizability and stability. We find a correlation between the Li ratio and the relative stability of CS versus non-CS frameworks, where the Li ratio is defined as the number of Li atoms divided by the number of non-Li metal atoms. Systems that prefer CS frameworks over non-CS frameworks generally exhibit small Li ratios (≤1.0). Fig. 3b illustrates this trend more clearly, showing a positive correlation between the Li ratio and E_hull^DFT (CS)–E_hull^DFT (non-CS). Even for several compositions with Li ratios greater than 1.0, no stable CS-topology structures are found (× markers in Fig. 3a).

Fig. 3 Structural stability of oxide SSEs with CS bond topology in representative element sets. (a) Upper part: Li ratio of the searched materials. Lower part: E_hull^DFT of CS and non-CS frameworks in each composition. Compositions where no stable CS structure is found are represented by × markers. The lowest-energy atomistic structures of LiTa₃P₂O₁₃, Li₂ZrSiO₅, Li₂AlPO₅, Li₆Al₂Si₃O₁₂, and Li₅Ti₂AlO₈ are presented below the bar plot as example structures at each Li ratio. (b) E_hull^DFT (CS)–E_hull^DFT (non-CS) as a function of Li ratio. (c) Structures of CS and non-CS frameworks in LiTaSi₂O₇ and Li₃AlSi₂O₇ and their energy above hull values (unit: meV atom⁻¹). (d) Ionic conductivity at 1000 K calculated by SevenNet-0 as a function of Li ratio for both CS and non-CS framework structures.

One might question our conclusion regarding the higher preference for CS frameworks in systems with low Li contents, since TOPIC explores CS-topology configurations more exhaustively than non-CS configurations in SPINNER, which is limited to searching unit cells with fewer atoms than those accessible in TOPIC. However, the same conclusion is reached when comparing CS and non-CS structures generated solely by SPINNER, as shown in Fig. S6 and S7, although some of the lowest-energy CS structures differ from those found in TOPIC (usually, TOPIC identifies more stable CS configurations than SPINNER). This supports that the preference for CS bond topology indeed exists in systems with low Li contents even though the existence of lower-energy non-CS polymorphs cannot be completely ruled out. Beyond our justification based on the simulation results presented above, the stronger preference for CS frameworks over non-CS ones at low Li ratios can be physically rationalized; crystal structures with fully connected CS frameworks exhibit relatively compact atomic packing, which imposes spatial constraints that limit the accommodation of Li ions. On the other hand, structures with high Li contents in well-known sulfide and chloride SSEs typically include isolated polyhedra (e.g., PS₄), providing greater internal space for Li accommodation.^59–61 This is exemplified by comparing LiTaSi₂O₇ (low Li ratio) and Li₃AlSi₂O₇ (high Li ratio). As shown in Fig. 3c, both structures share the same CS framework. For Li₃AlSi₂O₇, however, the CS polymorph is higher in hull energy than its non-CS counterpart, while for LiTaSi₂O₇, the CS polymorph remains more stable. Additionally, in Fig. 3b, frameworks incorporating larger Zr ions (86 pm) consistently exhibit lower E_hull^DFT (CS)–E_hull^DFT (non-CS) values than frameworks containing smaller Ti ions (75 pm). This trend can be attributed to the greater internal free space in the former, arising from the larger ionic radius of Zr⁴⁺ compared to Ti⁴⁺ (further illustrated in Fig. S8). The weak Li–O COHP, as discussed above, indicates low covalency and predominantly ionic interactions between Li and O ions, suggesting that Li influences structural stability mainly through electrostatic effects. Consequently, at low Li contents, the CS framework—composed of networks of directional covalent bonds between non-Li cations and oxygen—can accommodate Li ions with minimal structural distortion. As the Li content increases, however, enhanced Li⁺–Li⁺ electrostatic repulsion and the reduced availability of vacant sites for Li ions destabilize the compact corner-sharing framework, thereby explaining the observed dependence of structural stability on the Li ratio.

We examine the ionic conductivities of all lowest-energy and metastable structures within 50 meV atom⁻¹ above the hull at 1000 K. To this end, we perform MLIP-MD simulations using SevenNet-0, which has been shown to predict Li-ion conductivities with reasonable accuracy (see Fig. S9 and Table S3 for details).⁶² Fig. 3d shows the conductivity as a function of the Li ratio. Among 339 CS-framework structures examined, 133 (39%) exhibit conductivities above 101 mS cm⁻¹, which corresponds to practically relevant values of 0.1 mS cm⁻¹ at room temperature, assuming an activation energy of 0.3 eV.^1,43 In contrast, only 2 out of 70 non-CS structures (3%) exceed this threshold. As expected, frameworks with isolated vertices exhibit stronger Li–O interactions and hence lower ionic conductivity—consistent with the generally poorer transport of oxide frameworks with isolated vertices (e.g., oxide LISICON-type) relative to CS frameworks (e.g., NASICON-type).⁶³ Overall, low Li concentration favors the formation of CS-framework structures, which exhibit significantly higher Li-ion conductivity compared to non-CS structures. Therefore, we suggest that the Li ratio can serve as a key descriptor for screening oxide SSE materials.

2.5 Comprehensive screening of quaternary compositions favoring CS frameworks

We conduct a thorough screening of oxide SSE candidates in the quaternary compositional space using the Li-ratio descriptor. We first enumerate compositions of the form Li_x(OcO_6/2)_m(TO_4/2)_n, considering Oc : T ratios of 2 : 1, 3 : 2, 1 : 1, 2 : 3, and 1 : 2, under the constraint that the Li ratio is less than or equal to 1.0. The Oc and T elements discussed in the previous section are included in this screening. Furthermore, Mg, Ga, Sc, Hf, and Nb, which fulfill Pauling's rule but have been relatively underexplored, are included as Oc elements. This procedure yields a total of 45 quaternary compositions. For each composition, we apply both TOPIC and SPINNER and retain the lower-DFT-energy structure as described above. We observe a similar correlation between the Li ratio and E_hull^DFT (CS)–E_hull^DFT (non-CS) (Fig. S10 and Table S4) as in the previous screening results (Fig. 3b), indicating that the Li ratio descriptor is robust for these systems as well.

To efficiently identify candidates with high ionic conductivity, we employ a stepwise screening procedure, as illustrated in Fig. 4. First, we identify 45 quaternary Li-oxide compositions with Li ratios below 1.0, as described above. Next, we perform structure predictions with SPINNER (Z = 2 or 4, yielding unit cells with approximately 50 atoms, 60 000 structure generations) and TOPIC (Z = 4, 6, 300 000 trials for each). We retain only those with E_hull^DFT (CS) < E_hull^DFT (non-CS), reducing the set from 45 to 30. Subsequently, CS-topology crystal structures are generated with TOPIC from an expanded pool of 1.2 million trials with Z = 2, 4, 6, and 8 (300 000 trials for each), enabling more reliable identification of energetically favorable configurations. After discarding structures lying more than 100 meV atom⁻¹ above the lowest MLIP-energy for each composition and removing duplicates, we obtain 24 463 unique CS-framework structures. Finally, selecting those with E_hull^DFT ≤ 50 meV atom⁻¹ yields 438 distinct structures across 30 compositions. To rapidly evaluate ionic conductivity, a single MD simulation is performed for each selected configuration using SevenNet-0 at 1000 K, and only those with conductivities above 101 mS cm⁻¹ are retained, following ref. 43. For these structures, we conduct comprehensive MD simulations with 3–5 independent runs at each temperature (800, 900, 1000, 1100, and 1200 K). Ionic conductivities at 300 K are extrapolated using the Arrhenius relationship, as described in ref. 64. This screening identifies 92 candidates with extrapolated room-temperature ionic conductivities exceeding 0.1 mS cm⁻¹. Finally, for further validation, DFT-based MD simulations are carried out on 19 structures within 10 meV atom⁻¹ of the lowest energy for each composition (ΔE^DFT), ultimately yielding 15 candidate structures with room-temperature conductivities greater than 0.1 mS cm⁻¹. This energy range, ΔE^DFT, is accounted for to include potential candidates that are likely synthesizable due to their small energy differences, while simultaneously exhibiting high ionic conductivity, as demonstrated experimentally for LiAlSiO₄.^65,66 For the final candidates, we confirm that the conductivity values are not significantly affected by cell-size dependence (Fig. S11).

Fig. 4 Screening process. Orange boxes indicate target composition selection. Unique frameworks from these compositions are then evaluated as candidate solid electrolytes under several screening conditions (blue boxes).

Fig. 5 presents the 15 identified candidates, with detailed properties provided in Table 2 (see the Methods section for details of the calculations). Several of these candidates share frameworks with known compounds, including Li₂HfSiO₅(I) (LiTaSiO₅-type), Li₂HfSiO₅(II) (β-LiTiPO₅-type), LiNb₂PO₈(III) (LiTa₂PO₈-type), and three NASICON-type structures: Li₄Hf₂Si₃O₁₂(I), Li₄Zr₂Si₃O₁₂(I), and Li₂Nb₂Si₃O₁₂(II). Note that the Roman numerals in parentheses indicate the stability order among configurations of the same composition, corresponding to their energy ranking relative to the ground state (with I denoting the ground state). In addition, TOPIC discovers materials with novel frameworks. For example, TOPIC reveals a distinct structural feature, commonly observed in Li₂HfSi₂O₇, Li₂ZrSi₂O₇, and Li₃ScSi₂O₇, consisting of Si₂O₇ polyanions interconnected by OcO₆ octahedra. The conduction pathways in these materials are found to be quasi two-dimensional (see Fig. S12a). Another structural type, featuring similar Si₂O₇ and OcO₆ connectivity but a slightly different atomic arrangement, is shared by Li₂HfSi₂O₇ and Li₂ZrSi₂O₇ (see the blue dashed box in Fig. 5), leading to the three-dimensional conduction pathways (see Fig. S12b). LiNb₂PO₈(II) adopts a framework identical to that of metastable LiTa₂PO₈, newly identified by TOPIC in the previous subsection (Fig. 2b). Li conduction pathways of three additional structures with previously unreported frameworks—Li₂Nb₂Si₃O₁₂(I), LiNbSi₂O₇(VIII), and Li₃Nb₃Si₂O₁₃(X)—are also shown in Fig. S12. We find that several candidates in this set exhibit high ionic conductivities above 1 mS cm⁻¹. In particular, the NASICON-type Li₄Hf₂Si₃O₁₂ shows an ionic conductivity of 14.09 mS cm⁻¹ (7.56–26.22 mS cm⁻¹, considering uncertainty), which is comparable to that of liquid electrolytes and to the highest values reported for SSEs, such as LiNbOCl₄ and LiTaOCl₄.⁵⁸

Fig. 5 Final candidates. Structures are selected based on the criteria of high ionic conductivity (≥0.1 mS cm⁻¹) and a low energy difference (≤10 meV atom⁻¹) from the lowest-energy structure at each composition. The corresponding properties are listed in Table 2.

Table 2 Properties of candidate structures found by TOPIC. Energy above hull E_hull^DFT, energy difference from the lowest-energy structure at the given composition(ΔE^DFT), bandgap (E_g), activation energy (E_a), and ionic conductivity at 300 K (σ_300K) are provided. For E_a and σ_300K, estimates of the upper and lower bounds are shown in parentheses

Formula	Ehull^DFT (meV atom^−1)	ΔE^DFT (meV atom^−1)	Eg (eV)	Ea (eV)	σ300K (mS cm^−1)
Li2HfSi2O7 (I)	12	0	6.8	0.319 ± 0.070	0.392 (0.026, 5.879)
Li2HfSi2O7 (III)	18	6.2	6.9	0.234 ± 0.082	2.388 (0.100, 56.801)
Li2HfSiO5 (I)	19	0	7.0	0.299 ± 0.020	1.127 (0.516, 2.463)
Li2HfSiO5 (II)	19	0.2	6.9	0.375 ± 0.063	0.113 (0.010, 1.276)
Li4Hf2Si3O12 (I)	21	0	6.5	0.183 ± 0.016	14.090 (7.569, 26.229)
LiNb2PO8 (II)	25	9.3	4.2	0.353 ± 0.080	0.092 (0.004, 1.990)
LiNb2PO8 (III)	26	10.0	3.8	0.201 ± 0.034	7.820 (2.079, 29.416)
Li2ZrSi2O7 (VI)	28	2.1	6.3	0.277 ± 0.017	1.245 (0.657, 2.356)
Li2ZrSi2O7 (VII)	28	2.4	6.5	0.288 ± 0.067	0.553 (0.042, 7.335)
Li3ScSi2O7 (II)	28	6.0	6.7	0.236 ± 0.075	3.375 (0.188, 60.750)
LiNbSi2O7 (VIII)	29	9.2	5.0	0.328 ± 0.019	0.161 (0.078, 0.332)
Li2Nb2Si3O12 (I)	31	0	4.5	0.301 ± 0.085	0.232 (0.009, 6.211)
Li2Nb2Si3O12 (II)	32	0.6	4.8	0.267 ± 0.023	3.006 (1.245, 7.256)
Li4Zr2Si3O12 (I)	34	0	6.0	0.270 ± 0.006	0.756 (0.599, 0.955)
Li3Nb3Si2O13 (X)	44	9.7	3.9	0.274 ± 0.021	2.182 (0.959, 4.961)

---

We further evaluate the electrochemical stability of the final candidate materials by calculating their electrochemical stability windows (ESWs) following ref. 72 (see Fig. S13). Overall, the predicted candidates exhibit ESW ranges comparable to those of representative known oxide SSEs and generally wider than those of typical sulfide SSEs. We also note that the practically accessible stability range may be broader when indirect decomposition pathways are taken into account, as passivating interphases can partially mitigate interfacial reactions with electrodes.⁷³ Although this analysis provides only a preliminary filter, the present study does not include practical assessments of stability factors such as moisture tolerance and phase stability at finite temperature, which will need to be addressed in future studies.

In addition to the final candidates, we find that several structural frameworks (e.g., NASICON and new prototypes 1 and 3) recur across new materials generated by TOPIC that satisfy E_hull^DFT ≤ 50 meV, albeit with varying Li ratios (see Fig. 6a). Comparing the ionic conductivities of compounds sharing identical frameworks but different Li contents can provide crucial insights for discovering novel SSEs. In particular, stuffing Li into SSEs is a well-recognized strategy for enhancing Li-ion conductivity, as it lowers the migration barrier through increased Li–Li repulsion and the resulting structural distortions.^43,53 To examine whether this effect is also present in our candidates, we compare the activation barriers as a function of Li content for materials sharing the same framework. We analyze the NASICON-type, β-LiTiPO₅-type, and new prototypes 1 and 3, as shown in Fig. 6a. In all frameworks except new prototype 3, the activation barriers decrease as the Li ratio increases, which is consistent with previous theory. Indeed, when the activation barriers are plotted against the maximum CSM value for each prototype family, the barriers decrease as the CSM increases (Fig. S14). For instance, the activation barrier of synthesized β-LiTiPO₅ exceeds 1 eV, as measured experimentally,⁶⁷ whereas its analogue Li₂HfSiO₅(II), which has a higher Li ratio owing to the smaller valence of Si compared to P, exhibits a much lower barrier of 0.38 eV. CSM analysis supports this trend: Li sites in β-LiTiPO₅ show modest distortion (CSM = 2.2), while those in Li₂HfSiO₅(II) show significantly larger distortions (CSM = 5.4 and 8.2). In the atomistic structures in Fig. 6b, we also see that the Li sites are symmetric in β-LiTiPO₅, whereas they are distorted in Li₂HfSiO₅(II). The same trend is also observed in the NASICON-type and new prototype 1 structures (see Table S6 and Fig. S14).

Fig. 6 Effect of the Li ratio on ionic conductivity. (a) Activation energy as a function of Li ratio for NASICON-type, β-LiTiPO₅-type, new prototype 1, and new prototype 3. Markers indicate the average values of activation energy with the same framework and Li ratio, and error bars indicate the standard deviations. Activation energy values calculated with DFT are shown as filled markers, whereas experimentally reported values are shown as empty markers (empty diamond: β-LiTiPO₅-type;⁶⁷ empty circle: NASICON-type^68–71). The corresponding compounds for the NASICON-type and new prototype 3 structures are indicated in the figure, and details for the other types are provided in Table S5. (b) Structure and local environments of Li ions in β-LiTiPO₅ and Li₂TiSiO₅. CSM values and the site type (Oc for the octahedron and T for the tetrahedron) for each Li environment are provided below. (c) Li conduction pathways at 800 K for LiTa₂PO₈, Li₂Ta₂SiO₈, Hf-doped Li₂Ta₂SiO₈, and Li₃Ta₂AlO₈.

However, structures belonging to prototype 3 do not follow this trend. LiTa₂PO₈, despite having the lowest Li ratio, exhibits the lowest activation barrier because it contains quasi-one-dimensional channels that provide efficient conduction pathways. By contrast, in Li₂Ta₂SiO₈ (Li ratio = 0.667), the one-dimensional channels become overcrowded with Li ions, blocking not only the stable sites but also the transition-state sites along the migration pathway, which suppresses conductivity. Increasing the Li ratio further, as in Li₃Ta₂AlO₈ (Li ratio = 1.0), restores conductivity by opening alternative three-dimensional pathways. Similarly, introducing a small amount of aliovalent doping in Li₂Ta₂SiO₈—for example, substituting Ta with Hf—creates new three-dimensional pathways and thereby lowers the activation barriers from 0.562 eV to 0.352 eV. This trend can be explained by Li-site distortion: the activation barrier decreases as the CSM increases (Fig. S14), consistent with the tendency observed for the other prototype families. This result further supports the view that Li-site distortion is an important structural factor governing ionic conductivity.

To further examine which structural descriptors are related to ionic conductivity within diverse CS-framework structure prototypes reported in the previous study,⁴³ we compare polyhedra packing ratio, distant Li-site ratio, and CSM value between all 438 CS-framework structures generated by TOPIC with E_hull^DFT ≤ 50 meV atom⁻¹ and the 15 high-conductivity candidate structures (Fig. S15). However, these descriptors cannot capture conductivity trends across distinct prototype groups. Therefore, a more complete understanding of the relationship between structural distortion and ionic conductivity require more comprehensive future analysis.

3 Discussion

In this study, we have developed TOPIC for the prediction of CS-topology crystal structures in Li-oxide SSEs. For a single quaternary composition containing approximately 100 atoms per unit cell, the prediction of the lowest-energy structure requires 4 days of computational time. Given the intrinsic complexity of quaternary oxides and their large unit cells, this performance illustrates the efficiency of TOPIC, which requires only about twice the computational cost of earlier methods applied to simpler ternary oxides (about 50 atoms per unit cell).²⁴

By employing TOPIC, we have identified potential Li-oxide SSEs with high Li conductivity, as summarized in Table 2. We have discovered a diverse set of candidate materials exhibiting previously unreported structural prototypes and elemental combinations. Interestingly, while high-performance solid electrolytes containing silicon have been rarely reported in previous studies, our screening reveals that 13 out of 15 low-energy Si-containing compositions exhibit high predicted ionic conductivity. Among these, five materials contain Hf, which is also uncommon among reported compounds. Notably, Li₄Hf₂Si₃O₁₂ exhibits a high Li-ion conductivity of 14.09 mS cm⁻¹ (7.56–26.22 mS cm⁻¹, considering uncertainty), along with a wide band gap of 6.5 eV predicted using the accurate HSE06 method. In particular, this compound is predicted to be the lowest-energy structure among those with the same composition, implying relatively high synthesizability. Furthermore, we report various novel structural prototypes in this study (Fig. 5), which can serve as a materials library for designing new compounds through modifications such as phase distortions⁷⁴ and aliovalent doping.^75,76 For instance, we test the effectiveness of an aliovalent doping strategy on Li₂Ta₂SiO₈, whose lowest-energy structure with a novel framework exhibits low ionic conductivity (0.0004 mS cm⁻¹) at 300 K. With Hf doping (Li_2.125Ta_1.875Hf_0.125SiO₈), the activation energy decreases by 0.210 eV, and the ionic conductivity at 300 K increases by more than an order of magnitude (0.193 mS cm⁻¹). We also observe that doping induces a transition from a quasi-1D path in the undoped structure to a 3D Li conduction network, highlighting that doping creates additional conduction pathways (Fig. 6c). This is consistent with the transition to a 3D conduction pathway that occurs with increasing Li content across different element sets, as described in the previous section.

Even if this work conducts a comprehensive search of quaternary compositions composed of octahedral and tetrahedral polyhedra, it does not encompass the full chemical space of oxide SSEs. For example, the present screening does not include elements forming 3-coordinated polyhedra (e.g., SeO₃ in LiGa(SeO₃)₂ (ref. 43)) or garnet-type frameworks with La atoms such as LLZO. Furthermore, quinary compositions and doped materials may also exhibit high ionic conductivity. The structure libraries constructed in this work provide a valuable foundation for pursuing such extensions. The discovery of these novel compositions and prototypes is enabled by the ability of the TOPIC algorithm to directly explore the potential energy surface, thereby overcoming the limitations of conventional template-based methods.⁷⁷

Although the current implementation of TOPIC is specialized to oxide systems with fully connected CS frameworks, the framework-generation strategy could be extended. For instance, TOPIC can also be extended to mixed CS and isolated frameworks, corresponding to the broader definition of CS frameworks in ref. 43. In such cases, one could first generate a bond network with CS topology between neighboring cations and then place isolated anions around cations whose coordination numbers remain lower than their target values (see Fig. S16 for details). This type of extension may be particularly relevant for oxyhalides, where oxygen often participates in the connected framework while halogen atoms remain isolated. By contrast, extending TOPIC to structures containing both edge-sharing and corner-sharing motifs is expected to be more challenging, because the number of possible ways to assign sharing relationships between cations increases substantially. In fact, recent work has shown that, in sulfide SSEs, metastable edge-sharing frameworks can exhibit higher ionic conductivity than corner-sharing frameworks.³⁹ Thus, extending this type of analysis to broader chemistries and alternative framework connectivities will be an important direction for future work.

4 Conclusions

In summary, we develop a new structure prediction algorithm, TOPIC, designed to efficiently generate oxide SSE materials with CS frameworks. We validate that TOPIC successfully reproduces known oxide SSEs (up to 200 atoms per unit cell) as well as previously unreported polymorphs. By integrating a pretrained MLIP, SevenNet-0, we explore uncharted chemical spaces where the CS framework is likely to be preferred and identify several potential candidates. Our analysis reveals that the Li ratio plays a critical role in determining the stability of CS frameworks and provides a useful descriptor for screening. We analyze the effect of Li site distortion on ionic conductivity from the discovered structures, providing guidance for establishing design principles. Overall, TOPIC offers a scalable route to explore uncharted chemical space for next-generation oxide electrolytes.

5 Methods

5.1 TOPIC algorithm

5.1.1 Random structure generation for cation sublattices. Cation sublattices are generated under specified space groups using the RandSpg package,⁴⁶ which randomly assigns lattice vectors and Wyckoff positions while preserving the designated symmetry. The symmetry group of each cell is chosen randomly. The mean cell volumes of the generated structures are referenced to the density of the amorphous phase obtained from DFT melt–quench simulations (see below). Volumes are selected within 90–120% of the mean cell volume, with each lattice parameter constrained to 40–250% of the cubic root of the mean cell volume. Lattice angles are restricted to 60°–120°. Distance constraints are applied to all cation–cation pairs by imposing a minimum separation equal to the values observed for the corresponding pairs during the melt–quench–annealing trajectories. This procedure ensures chemically reasonable configurations and avoids unphysical short interatomic distances.

5.1.2 Oxygen placement. After generating random cation sites, oxygen atoms are placed at the midpoints of cation pairs. Two different algorithms are employed to generate cation pairings (virtual bonds) that satisfy both the CS bond topology and coordination number constraints. (1) In the first algorithm, the bond topology is generated by sequentially forming virtual bonds for each cation according to its predefined coordination number. Specifically, the order of cations is first randomly determined. For each cation in this order, connections are made to the nearest neighboring cations until its coordination number is satisfied. If a virtual bond has already been assigned by another cation, the bond is retained, and additional connections are made preferentially to the nearest unsaturated neighbors. If the nearest cation has already reached its coordination number limit, the next nearest cation is chosen, and this process is repeated until all cations satisfy their coordination numbers. When this process fails (e.g., by forming isolated loops), a new random cation order is selected, and the procedure is retried up to 1000 times. If all attempts fail, or if the generated structure succeeds but has P1 symmetry, a new random cation configuration is generated and the algorithm is repeated until success. Once the bond topology is completed, oxygen atoms are placed at the midpoints of the virtual bonds. (2) In the second algorithm, all cation pairs with short separations that can form virtual bonds are first enumerated, and only those structures in which a CS bond topology is generated are selected. Specifically, all cation pairs with separations shorter than twice the cation–O distance (the first radial distribution function (RDF) valley from the MD simulations) are considered, and oxygen atoms are placed at the midpoints accordingly. After assigning oxygen atoms, the coordination number of each cation is checked. If any cation does not satisfy its predefined coordination number, the cation structures are regenerated under the same space group, and this process is repeated up to 100 trials. If all attempts fail, a new space group number is randomly selected and the process is repeated until success. Compared with the first method, this approach generally preserves the symmetry of the initial cation arrangement, and the resulting structures therefore tend to exhibit higher symmetry. In practice, structures are generated such that 50% are obtained from the first method and 50% from the second.

5.1.3 Bond-length optimization with classical potentials. Oxygen atoms are initially placed at the midpoints of the virtual bonds to satisfy the CS bond topology. However, the resulting cation–O bond lengths are not optimal. Using the LAMMPS package,⁷⁸ the structures are optimized using a Lennard-Jones (LJ) potential,
which is selectively applied to the virtual bonds between cations and oxygen atoms. Note that the LJ potentials are not applied to all pairs of atoms, but are selectively applied to O–cation pairs connected by virtual bonds. Here, σ is set to the position of the first peak of the cation–O RDF obtained from amorphous structures generated by DFT-MD simulations for the training set, divided by , ensuring that the equilibrium bond length coincides with the RDF peak. The ε value is set to 3.0 eV for all cases. To prevent the formation of highly unphysical structures, additional harmonic repulsion potentials are introduced,

E = K(R − R_c)²

whenever two atoms are closer than the cutoff distance R_c (fix restrain command in LAMMPS). The K value is set to 0.2 eV Å⁻² for all cases. These repulsion terms are applied between O–O pairs, between cation–cation pairs, and also between cations and O atoms belonging to different polyhedra. For the cutoff distances of O–O and cation–cation pairs, the shortest distances observed in melt–quench–annealing MD trajectories are used. In the case of cation–oxygen repulsion across different polyhedra, the cutoff is set to the position of the first valley in the cation–O RDF of the amorphous structure.

5.1.4 Optimizing frameworks. After optimization with the LJ potential, the code checks whether the structure satisfies the desired coordination numbers and CS bond topology. If this condition is not met, the structure is discarded; otherwise, it is further relaxed with the MLIPs trained without labeling the Li atoms. The details of training set generation and model architectures are described below. The structure is first optimized using an MLIP with a cutoff of 3 Å, and then further relaxed using a model with a cutoff of 6.0 Å. At each stage, structures that do not preserve the CS bond topology or bond-length conditions after relaxation are discarded.

5.1.5 Li insertion process. We use Voronoi tessellation based on oxygen atom positions⁴³ to enumerate possible Li sites using the SciPy package.⁷⁹ Among the Voronoi nodes, those located within 1 Å of non-Li cations are discarded. The remaining nodes are clustered if they are within 1 Å of each other and represented by the average position. To obtain stable configurations, Li atoms are first placed in the largest free spaces and the structures are relaxed using MLIPs (trained for all elements including Li). If the structure collapses such that the coordination numbers of the non-Li cations are altered, Li atoms are inserted at random candidate sites and the structures are checked again. If no structure satisfying the conditions is obtained within 10 trials, the framework is considered unstable upon Li insertion and discarded. If the framework remains stable, 0 K Monte Carlo (MC) simulations are performed. In each MC step, one Li atom is randomly selected and its position is changed to one of the candidate Li sites. After relaxation with MLIPs, the bond topology and coordination numbers are checked. If the conditions are not satisfied or the energy increases, the move is rejected. After 100 MC steps, the lowest-energy structure is selected.

5.1.6 Final energy evaluations. After generating candidate structures with the above processes, the final candidates within an energy window of 100 meV atom⁻¹ from the MLIP calculations are selected and then recalculated using DFT. To improve efficiency, we adopt a two-step procedure. First, one-shot calculations are performed for all candidate structures using a k-point grid with the same k-spacing as in the DFT melt–quench simulations. Next, structures within an energy window of 50 meV atom⁻¹ from the one-shot DFT results are selected and relaxed with the AMP² package⁸⁰ to obtain accurate energies. Convergence criteria for the k-point tests are set to 3 meV atom⁻¹ for energy and 10 kbar for pressure. The cutoff energy for all DFT calculations is set to 520 eV.

5.1.7 Modification of the TOPIC algorithm for Li₇La₃Zr₂O₁₂. Generating the (ZrO_6/2)₂(LiO_4/2)₃ framework follows the standard TOPIC workflow: sampling random cation sublattices and optimizing bond lengths with LJ potentials. During framework optimization, however, we use only the 3 Å-cutoff MLIP because the 6 Å-cutoff model exhibits roughly twice the force error of our validation systems, whereas the 3 Å model is only 40% higher. During the La-insertion MC, we evaluate single-point energies rather than performing relaxations in order to avoid artifacts arising from the non-stoichiometric conditions. Finally, we remove the Li atoms within the framework and reinsert Li atoms via MC simulations, yielding the final Li₇La₃Zr₂O₁₂ structures.

5.2 Machine learning interatomic potential training

5.2.1 MLIP details. We employ Behler–Parrinello neural network potentials,³² with symmetry function vectors as descriptors, using the SIMPLE-NN package.⁸¹ For each atomic pair, 8 radial and 18 angular symmetry functions are constructed with identical cutoff radii applied to both. Model training is continued until the root-mean-square errors (RMSEs) of the validation set fail to decrease by at least 0.5 meV atom⁻¹ (energy), 0.01 eV Å⁻¹ (force), or 0.5 kbar (stress) over 50 epochs. The average validation RMSEs are: (1) 39 meV atom⁻¹ (energy), 2.0 eV Å⁻¹ (force), and 59 kbar (stress) for the 3 Å cutoff models excluding Li atoms; (2) 14 meV atom⁻¹, 0.87 eV Å⁻¹, and 44 kbar for the 6 Å cutoff models excluding Li atoms; and (3) 4 meV atom⁻¹, 0.33 eV Å⁻¹, and 8 kbar for the 6 Å cutoff models including all atoms. The neural networks consist of two hidden layers with 30 nodes each. Input vectors are decorrelated using principal component analysis (PCA) and subsequently whitened to accelerate convergence.⁸² Model weights are optimized using the Adam algorithm,⁸³ with a batch size of 20. To mitigate overfitting, L2 regularization is incorporated into the loss function, and 10% of the dataset is reserved for validation.

5.2.2 Generating training set. Disordered structures for the training set are generated using melt–quench–annealing MD simulations with DFT calculations. Initial configurations are constructed by randomly placing atoms within a volume estimated from their atomic radii. The system is first equilibrated at 4500 K for 3 ps, followed by melting at an empirically determined melting temperature for 8 ps. The melting temperature is estimated as the point at which the diffusion coefficient of all atomic species exceeds 4 × 10⁻⁹ m² s⁻¹, a common characteristic value for liquid metals and ceramics.^84,85 From this molten state, the system is quenched at a cooling rate of 200 K ps⁻¹ and subsequently annealed at 500 K for 4 ps, yielding amorphous configurations suitable for MLIP training and providing preliminary structural information.

5.2.3 Density functional theory calculations. All DFT calculations are performed using the Vienna Ab initio Simulation Package (VASP)⁸⁶ with projector augmented-wave (PAW) pseudopotentials⁸⁷ and the Perdew–Burke–Ernzerhof (PBE) functional⁸⁸ for the exchange–correlation energy. The cutoff energies are determined from convergence tests on premelted structures, ensuring that the total energies, forces, and stress tensors converge within 20 meV atom⁻¹, 0.3 eV Å⁻¹, and 10 kbar, respectively. MLIPs trained on melt–quench–anneal trajectories that satisfy these thresholds are shown to be sufficiently accurate for crystal-structure prediction, successfully generating low-energy structures in the process.^24,40 A single k-point, either Γ or ,⁸⁹ is used for the Brillouin-zone integration, with the choice determined by the same convergence tests as for the cutoff energies. Band gap values of the final candidates are evaluated by one-shot hybrid functional calculations (HSE06 functional)⁹⁰ on the PBE-optimized structures using the AMP² package.⁹¹

5.3 SPINNER

We use SPINNER to generate structures without bond topology constraints.⁴⁰ Structures are produced through random seeding, permutation, and lattice mutation in ratios of 70, 20, and 10%, respectively. The population size of each generation is capped at 300, and the search proceeds for 200 generations. In each generation, structures within a 100 meV atom⁻¹ window from the lowest-energy structure of the previous generation are selected for mutation operations. After the entire run, the final candidate structures within the lowest 50 meV atom⁻¹ energy window are retained. During structural relaxations, we use the same MLIPs as in TOPIC, trained with a 6 Å cutoff radius and including all atoms.

5.4 Li ionic conductivity calculation

We follow the method described in ref. 64. MD simulations up to 50 ps (both DFT and MLIP) are conducted 3–5 times, and a time-averaging method is used to calculate the mean squared displacement (MSD), defined as follows:where N_Li is the number of Li atoms, R_i(t) is the position vector of Li atom i at time t, and τ is the lag time. The notation 〈⋯ 〉_t denotes an average over all permissible time origins t within a trajectory. The diffusion coefficient at each temperature is calculated from the Einstein relation:where D₀ is the pre-factor of the diffusion coefficient, E_a is the activation energy of diffusion, k_B is the Boltzmann constant, and T is the temperature. To obtain a linear region of MSD, we remove the initial 2 ps of MSD and beyond 50% of the total simulation time. We verify that all MSD curves obtained from DFT MD of the 15 candidate structures in this study are within the linear region. The activation energy and ionic conductivity are obtained by fitting the diffusion coefficients to the Arrhenius equation in the temperature range of 800–1200 K:where C is a temperature-independent constant. We apply a weighted linear fit to the Arrhenius relation using weights of the inverse square of the standard deviation of log(D_Li) to each data point to handle the statistical uncertainty and thereby estimate the activation energy and its uncertainty.⁹²

The ionic conductivity is then calculated from the Nernst–Einstein relation:

where e is the elementary charge and V is the volume of the simulated cell. We set z = 1 for Li ions. The uncertainty in the ionic conductivity was estimated by propagating the uncertainty in the activation energy through the Arrhenius relation.

For visualization of Li conduction pathways, the pymatgen-analysis-diffusion package is employed,^93,94 and plotted isosurfaces of Li probability density at a threshold of 0.001 P_max, where P_max denotes the maximum value of the Li-ion probability density distribution.

While SevenNet-0 shows good predictive performance for materials with known structures,^95–97 its reliability for newly discovered prototypes remains uncertain. To assess its applicability to evaluate Li-ion transport, we compare diffusion coefficients predicted by SevenNet-0 with those obtained from DFT (Fig. S17). Overall, SevenNet-0 does not show the clear uniform shift relative to DFT that would be expected for a severe systematic softening of the potential energy surface (Fig. S17a).⁶² Instead, the deviations are strongly structure-dependent. For compounds with known prototypes, the model achieves a mean percentage error (MPE) of 8% and a mean absolute percentage error (MAPE) of 23%. In contrast, for the novel frameworks discovered in this study, SevenNet-0 tends to overestimate the diffusion coefficients, with an MPE of 26% and a MAPE of 38%, likely arising from a softening of the potential energy surface relative to DFT.^62,98 Nevertheless, these errors do not alter the order of magnitude of the predicted conductivities. Thus, SevenNet-0 remains suitable for high-throughput screening of oxide-based Li-ion solid electrolytes.

5.5 Aliovalent doping

Voronoi nodes are considered as candidate Li sites, as described above. For each system, 50 structures with random cation substitutions and Li stuffing are generated. Their energies are evaluated using SevenNet-0, and the lowest-energy structure is further relaxed with SevenNet-0. Although the selected structures may not correspond to the absolute lowest-energy configurations, they are expected to represent the doped systems reasonably well, as the ionic conduction characteristics are reproduced in line with previous reports (e.g., Ga-doped LiTiPO₅ and P-doped Li₂Mg₂(SO₄)₃, see Table S3).

5.6 Framework novelty assessment

To assess the novelty of the crystal structures identified in this work, we compare them against the AFLOW prototype library.⁹⁹ For each candidate, we examine both the full structure including Li atoms and the corresponding framework excluding Li. If neither structure matches any entry in the AFLOW prototype library, we classify the framework as novel.

Author contributions

S. K. conceived the original idea of topology constraints and developed the prototype code. S. Hw. developed the main algorithm of TOPIC and carried out the major tasks, including code implementation, MLIP training, performance tests, screening, and data analysis. J. L. trained MLIPs for a subset of compositions. J. K. contributed to analyzing diffusion properties. S. K., Y. K., and S. Ha. supervised the project and provided discussions. All authors contributed to drafting the manuscript.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The main part of TOPIC is opened at (https://github.com/kang1717/TOPIC). SPINNER for carrying evolutionary structure searches is available at (https://github.com/MDIL-SNU/SPINNER).⁴⁰ SIMPLE-NN for training MLIPs is available at (https://github.com/MDIL-SNU/SIMPLE-NN_v2).⁸¹ The SevenNet code utilized in this study is available in the project GitHub repository (https://www.github.com/MDIL-SNU/SevenNet).⁴⁵

The energy values, diffusion properties, and corresponding structure files for the discovered structures are provided separately as supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6ta01092j.

Acknowledgements

This research was supported by the Nano & Material Technology Development Programs through the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (RS-2024-00407995 and RS-2024-00450102). The computations were carried out at the Korea Institute of Science and Technology Information (KISTI) National Supercomputing Center (KSC-2025-CRE-0110).

References

1. Z. Zhang, Y. Shao, B. Lotsch, Y.-S. Hu, H. Li, J. Janek, L. F. Nazar, C.-W. Nan, J. Maier, M. Armand and L. Chen, Energy Environ. Sci., 2018, 11, 1945–1976.
2. A. C. Luntz, J. Voss and K. Reuter, J. Phys. Chem. Lett., 2015, 6, 4599–4604.
3. A. Manthiram, X. Yu and S. Wang, Nat. Rev. Mater., 2017, 2, 16103.
4. W. Li, J. R. Dahn and D. S. Wainwright, Science, 1994, 264, 1115–1118.
5. Q. Wang, P. Ping, X. Zhao, G. Chu, J. Sun and C. Chen, J. Power Sources, 2012, 208, 210–224.
6. J. Lau, R. H. DeBlock, D. M. Butts, D. S. Ashby, C. S. Choi and B. S. Dunn, Adv. Energy Mater., 2018, 8, 1800933.
7. R. Chen, Q. Li, X. Yu, L. Chen and H. Li, Chem. Rev., 2020, 120, 6820–6877.
8. K. J. Kim, M. Balaish, M. Wadaguchi, L. Kong and J. L. M. Rupp, Adv. Energy Mater., 2021, 11, 2002689.
9. R. Murugan, V. Thangadurai and W. Weppner, Angew. Chem., Int. Ed., 2007, 46, 7778–7781.
10. S. Stramare, V. Thangadurai and W. Weppner, Chem. Mater., 2003, 15, 3974–3990.
11. H. Aono, E. Sugimoto, Y. Sadaoka, N. Imanaka and G.-y. Adachi, J. Electrochem. Soc., 1990, 137, 1023–1027.
12. J. Kim, J. Kim, M. Avdeev, H. Yun and S.-J. Kim, J. Mater. Chem. A, 2018, 6, 22478–22482.
13. G. Bergerhoff, R. Hundt, R. Sievers and I. D. Brown, J. Chem. Inf. Comput. Sci., 1983, 23, 66–69.
14. X. He, Q. Bai, Y. Liu, A. M. Nolan, C. Ling and Y. Mo, Adv. Energy Mater., 2019, 9, 1902078.
15. L. Kahle, A. Marcolongo and N. Marzari, Energy Environ. Sci., 2020, 13, 928–948.
16. B. He, S. Chi, A. Ye, P. Mi, L. Zhang, B. Pu, Z. Zou, Y. Ran, Q. Zhao, D. Wang, W. Zhang, J. Zhao, S. Adams, M. Avdeev and S. Shi, Sci. Data, 2020, 7, 151.
17. A. D. Sendek, Q. Yang, E. D. Cubuk, K.-A. N. Duerloo, Y. Cui and E. J. Reed, Energy Environ. Sci., 2016, 10, 306–320.
18. K. Fujimura, A. Seko, Y. Koyama, A. Kuwabara, I. Kishida, K. Shitara, C. A. J. Fisher, H. Moriwake and I. Tanaka, Adv. Energy Mater., 2013, 3, 980–985.
19. Z. Zhu, I.-H. Chu and S. P. Ong, Chem. Mater., 2017, 29, 2474–2484.
20. J. Kim, D. H. Mok, H. Kim and S. Back, ACS Appl. Mater. Interfaces, 2023, 15, 52427–52435.
21. A. R. Oganov, C. J. Pickard, Q. Zhu and R. J. Needs, Nat. Rev. Mater., 2019, 4, 331–348.
22. C. W. Glass, A. R. Oganov and N. Hansen, Comput. Phys. Commun., 2006, 175, 713–720.
23. S. T. Call, D. Y. Zubarev and A. I. Boldyrev, J. Comput. Chem., 2007, 28, 1177–1186.
24. S. Hwang, J. Jung, C. Hong, W. Jeong, S. Kang and S. Han, J. Am. Chem. Soc., 2023, 145, 19378–19386.
25. A. R. Oganov, J. Chen, C. Gatti, Y. Ma, Y. Ma, C. W. Glass, Z. Liu, T. Yu, O. O. Kurakevych and V. L. Solozhenko, Nature, 2009, 460, 292.
26. D. Duan, Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao, H. Yu, B. Liu, W. Tian and T. Cui, Sci. Rep., 2014, 4, 6968.
27. W. Ming, M. Yoon, M.-H. Du, K. Lee and S. W. Kim, J. Am. Chem. Soc., 2016, 138, 15336–15344.
28. X. Wang, R. Xiao, H. Li and L. Chen, Phys. Rev. Lett., 2017, 118, 195901.
29. Y. Wang, Z. Ren, J. Zhang, S. Lu, C. Hua, H. Yuan, J. Luo, Y. Liu, J. Nai and X. Tao, Adv. Sci., 2024, 11, 2404213.
30. Q. Zhang, A. Choudhury and A. Chernatynskiy, J. Phys.: Condens. Matter, 2025, 37, 095901.
31. J. Clymo, C. M. Collins, K. Atkinson, M. S. Dyer, M. W. Gaultois, V. V. Gusev, M. J. Rosseinsky and S. Schewe, Angew. Chem., Int. Ed., 2025, 137, e202417657.
32. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 98, 146401.
33. A. V. Shapeev, Multiscale Model. Simul., 2016, 14, 1153–1173.
34. A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Phys. Rev. Lett., 2010, 104, 136403.
35. E. V. Podryabinkin, E. V. Tikhonov, A. V. Shapeev and A. R. Oganov, Phys. Rev. B, 2019, 99, 064114.
36. C. Hong, J. M. Choi, W. Jeong, S. Kang, S. Ju, K. Lee, J. Jung, Y. Youn and S. Han, Phys. Rev. B, 2020, 102, 224104.
37. Q. Tong, L. Xue, J. Lv, Y. Wang and Y. Ma, Faraday Discuss., 2018, 211, 31–43.
38. V. L. Deringer, C. J. Pickard and G. Csányi, Phys. Rev. Lett., 2018, 120, 156001.
39. J. H. Kim, J. S. Kim, Y. H. Kim, B. Jun, Y. J. Jang and S. U. Lee, J. Am. Chem. Soc., 2025, 147, 47381–47391.
40. S. Kang, W. Jeong, C. Hong, S. Hwang, Y. Yoon and S. Han, npj Comput. Mater., 2022, 8, 108.
41. T. Sato, T. Otani, S. Nakamori, F. Utsuno, T. Honma and T. Yamashita, Jpn. J. Appl. Phys., 2024, 63, 075503.
42. Y. Han, C. Ding, J. Wang, H. Gao, J. Shi, S. Yu, Q. Jia, S. Pan and J. Sun, Nat. Comput. Sci., 2025, 5, 255–267.
43. K. Jun, Y. Sun, Y. Xiao, Y. Zeng, R. Kim, H. Kim, L. J. Miara, D. Im, Y. Wang and G. Ceder, Nat. Mater., 2022, 21, 924–931.
44. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
45. Y. Park, J. Kim, S. Hwang and S. Han, J. Chem. Theory Comput., 2024, 20, 4857–4868.
46. P. Avery and E. Zurek, Comput. Phys. Commun., 2017, 213, 208–216.
47. C. J. Pickard and R. J. Needs, Phys. Rev. Lett., 2006, 97, 045504.
48. C. J. Pickard and R. J. Needs, J. Phys.: Condens. Matter, 2011, 23, 053201.
49. L. Pauling, J. Am. Chem. Soc., 1929, 51, 1010–1026.
50. R. D. Shannon, Acta Crystallogr., Sect. A, 1976, 32, 751–767.
51. T. Böger, T. Bernges, M. T. Agne, P. Canepa, F. Tietz and W. G. Zeier, J. Am. Chem. Soc., 2024, 146, 32678–32688.
52. V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Phys. Chem. A, 2011, 115, 5461–5466.
53. X. He, Y. Zhu and Y. Mo, Nat. Commun., 2017, 8, 15893.
54. T. Hupfer, E. Bucharsky, K. Schell and M. Hoffmann, Solid State Ionics, 2017, 302, 49–53.
55. P. Norby, Zeolites, 1990, 10, 193–199.
56. M. Umair, S. Zhou, W. Li, H. T. H. Rana, J. Yang, L. Cheng, M. Li, S. Yu and J. Wei, Batter. Supercaps, 2024, 8, e202400667.
57. F. Hussain, P. Li and Z. Li, J. Phys. Chem. C, 2019, 123, 19282–19287.
58. Y. Tanaka, K. Ueno, K. Mizuno, K. Takeuchi, T. Asano and A. Sakai, Angew. Chem., Int. Ed., 2023, 62, e202217581.
59. H.-J. Deiseroth, S.-T. Kong, H. Eckert, J. Vannahme, C. Reiner, T. Zaiß and M. Schlosser, Angew. Chem., 2008, 120, 767–770.
60. N. Kamaya, K. Homma, Y. Yamakawa, M. Hirayama, R. Kanno, M. Yonemura, T. Kamiyama, Y. Kato, S. Hama, K. Kawamoto and A. Mitsui, Nat. Mater., 2011, 10, 682–686.
61. T. Asano, A. Sakai, S. Ouchi, M. Sakaida, A. Miyazaki and S. Hasegawa, Adv. Mater., 2018, 30, 1803075.
62. J. Kim, J. Lee, S. Oh, Y. Park, S. Hwang, S. Han, S. Kang and Y. Kang, npj Comput. Mater., 2025, 12, 26.
63. B. Tao, C. Ren, H. Li, B. Liu, X. Jia, X. Dong, S. Zhang and H. Chang, Adv. Funct. Mater., 2022, 32, 2203551.
64. X. He, Y. Zhu, A. Epstein and Y. Mo, npj Comput. Mater., 2018, 4, 18.
65. F. Shin-ichi, S. Satoshi, S. Kaduhiro and T. Hitoshi, Solid State Ionics, 2004, 167, 325–329.
66. R. Sheil, Y.-C. Perng, J. Mars, J. Cho, B. Dunn, M. F. Toney and J. P. Chang, ACS Appl. Mater. Interfaces, 2020, 12, 56935–56942.
67. A. Robertson, J. G. Fletcher, J. M. S. Skakle and A. R. West, J. Solid State Chem., 1994, 109, 53–59.
68. H. Aono, E. Sugimoto, Y. Sadaoka, N. Imanaka and G.-y. Adachi, Solid State Ionics, 1991, 47, 257–264.
69. K. Takada, M. Tansho, I. Yanase, T. Inada, A. Kajiyama, M. Kouguchi, S. Kondo and M. Watanabe, Solid State Ionics, 2001, 139, 241–247.
70. K. Arbi, M. A. París and J. Sanz, J. Phys. Chem. B, 2006, 110, 6454–6457.
71. Y. Luo, X. Liu, C. Wen, T. Ning, X. Jiang and A. Lu, Appl. Phys. A, 2023, 129, 518.
72. Y. Zhu, X. He and Y. Mo, ACS Appl. Mater. Interfaces, 2015, 7, 23685–23693.
73. T. K. Schwietert, V. A. Arszelewska, C. Wang, C. Yu, A. Vasileiadis, N. J. J. D. Klerk, J. Hageman, T. Hupfer, I. Kerkamm, Y. Xu, E. V. D. Maas, E. M. Kelder, S. Ganapathy and M. Wagemaker, Nat. Mater., 2020, 19, 428–435.
74. J. Awaka, N. Kijima, H. Hayakawa and J. Akimoto, J. Solid State Chem., 2009, 182, 2046–2052.
75. G. T. Hitz, E. D. Wachsman and V. Thangadurai, J. Electrochem. Soc., 2013, 160, A1248–A1255.
76. S. Adams and R. P. Rao, J. Mater. Chem., 2011, 22, 1426–1434.
77. D. W. Davies, K. T. Butler, J. M. Skelton, C. Xie, A. R. Oganov and A. Walsh, Chem. Sci., 2017, 9, 1022–1030.
78. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. i. Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
79. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. v. d. Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. v. Mulbregt, A. Vijaykumar, A. P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström, C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavič, J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. V. D. M. Cardoso, J. Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y. O. Halchenko and Y. Vázquez-Baeza, Nat. Methods, 2020, 17, 261–272.
80. Y. Youn, M. Lee, C. Hong, D. Kim, S. Kim, J. Jung, K. Yim and S. Han, Comput. Phys. Commun., 2020, 256, 107450.
81. K. Lee, D. Yoo, W. Jeong and S. Han, Comput. Phys. Commun., 2019, 242, 95–103.
82. D. Yoo, K. Lee, W. Jeong, D. Lee, S. Watanabe and S. Han, Phys. Rev. Mater., 2019, 3, 093802.
83. D. P. Kingma and J. Ba, arXiv, 2014, preprint, arXiv:1412.6980,  DOI:10.48550/arXiv.1412.6980.
84. P. Protopapas, H. C. Andersen and N. A. D. Parlee, J. Chem. Phys., 1973, 59, 15–25.
85. J. W. Nowok, J. Mater. Sci., 1996, 31, 5169–5177.
86. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
87. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
88. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
89. A. Baldereschi, Phys. Rev. B, 1972, 7, 5212–5215.
90. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
91. S. Kim, M. Lee, C. Hong, Y. Yoon, H. An, D. Lee, W. Jeong, D. Yoo, Y. Kang, Y. Youn and S. Han, Sci. Data, 2020, 7, 387.
92. D. Ruppert and M. P. Wand, Ann. Math. Stat., 1994, 22, 1346–1370.
93. Z. Deng, Z. Zhu, I.-H. Chu and S. P. Ong, Chem. Mater., 2017, 29, 281–288.
94. Z. Zhu, I.-H. Chu, Z. Deng and S. P. Ong, Chem. Mater., 2015, 27, 8318–8325.
95. A. Loew, D. Sun, H.-C. Wang, S. Botti and M. A. L. Marques, npj Comput. Mater., 2025, 11, 178.
96. J. Riebesell, R. E. A. Goodall, P. Benner, Y. Chiang, B. Deng, G. Ceder, M. Asta, A. A. Lee, A. Jain and K. A. Persson, Nat. Mach. Intell., 2025, 7, 836–847.
97. B. Póta, P. Ahlawat, G. Csányi and M. Simoncelli, arXiv, 2025, preprint, arXiv:2408.00755,  DOI:10.48550/arXiv.2408.00755.
98. B. Deng, Y. Choi, P. Zhong, J. Riebesell, S. Anand, Z. Li, K. Jun, K. A. Persson and G. Ceder, npj Comput. Mater., 2025, 11, 9.
99. D. Hicks, C. Toher, D. C. Ford, F. Rose, C. D. Santo, O. Levy, M. J. Mehl and S. Curtarolo, npj Comput. Mater., 2021, 7, 30.

---