Accounting for the vibrational contribution to the configurational entropy in disordered solids with machine learned forcefields: a case study of garnet electrolyte Li₇La₃Zr₂O₁₂

Abstract

Accounting for lattice vibrations to accurately determine the phase stabilities of site-disordered solids is a long-standing challenge in computational material designs, due to the high computational cost associated with sampling the vast configurational space to obtain the converged thermodynamic quantities. One example is the garnet electrolyte Li₇La₃Zr₂O₁₂, the high-temperature and high-ion-mobility cubic phase of which is disordered in its Li⁺ site occupations, such that both the vibrational and configurational entropic contributions to its phase stability cannot be ignored. Understanding the subtle interplay between vibrational and configurational entropies in this material will therefore play a critical role in the rational manipulation of dopants and defects to stabilise cubic Li₇La₃Zr₂O₁₂ at room temperature for practical applications. Here, by developing machine learned forcefields based on an equivariant message-passing neural network SO3KRATES, we follow a strict statistical thermodynamic protocol to quantify the phase stability of cubic Li₇La₃Zr₂O₁₂ through structural optimisations, as well as molecular dynamics simulations at 300 and 1500 K, for a total of 70 120 configurations of cubic Li₇La₃Zr₂O₁₂. Although this only covers a tiny fraction of the configurational space (∼7 × 10³⁴ configurations in total), we are able to deterministically show that the vibrational contributions to the total configurational free energy at 1500 K are significant (on the order of 1 eV per atom) in correctly ordering the stability of the cubic Li₇La₃Zr₂O₁₂ over its tetragonal counterpart, thanks to the high data efficiency, accuracy, stability and good transferability of the transformer-based equivariant network architecture behind SO3KRATES. Therefore, our work opens up new avenues to accelerate the accurate computational designs of disordered solids, such as solid electrolytes, for technologically important applications.

---

1 Introduction

Computational modelling is an essential component in driving modern scientific discoveries and innovations across multiple fields of physics, chemistry, biology and materials science. These models provide us with the critical insights into the atomistic, and often electronic, origins behind the exotic properties of materials, such as thermal,^1,2 optical,^3,4 catalytic⁵ properties and many others. Regardless of which model is employed, fore-mostly, it determines the potential energy of a material system as a function of the atomic positions E(R) with a given composition, a key quantity useful for determining a material's thermodynamic stability, thus providing an answer to the important question of whether the material can exist in the first place. The capability of these energy models, particularly density functional theory (DFT), to reach sub-kcal mol⁻¹ (in chemistry) or meV per atom (in physics and materials science) accuracies compared to experimental thermochemical data as demonstrated by many carefully curated benchmarks^6–10 has brought confidence in applying these models as the predictive tools to guide experimental material discoveries. A notable example is the application of organic^11,12 and inorganic^13,14 crystal structure predictions, which sample the potential energy surface (PES) for structures corresponding to the global energy minimum, and/or a set of low energy local minima, with demonstrated success in materials discoveries such as highly porous organic solids for gas adsorptions,¹⁵ complex layered inorganic solids,¹⁶ and new ternary inorganic nitrides.¹⁷

Despite these successes, structure predictions remain one of the Holy Grails in computational chemistry and materials science. This is simply because the size of the search space grows combinatorially large as the degrees-of-freedoms increase, making it impossible to conduct an exhaustive search over the PES, risking important structures being missed¹⁸ and/or difficulty in converging thermodynamic properties through ensemble averaging.¹⁹ This challenge is often compounded by the necessity to include the (vibrational) free energy contributions for accurate modelling of materials' thermodynamic stability at finite temperatures.^7,20–22

Garnet Li₇La₃Zr₂O₁₂ (LLZO) is an interesting example, and is the focus of this work, for demonstrating the challenges behind sampling the PES using computational approaches for comprehensively understanding the phase behaviours of disordered solids at the atomistic level [Fig. 1]. LLZO is a promising solid electrolyte for all solid-state lithium ion batteries because of its high Li⁺ ion conductivity²³ with negligible electronic transports,²⁴ wide electrochemical windows,²⁵ and chemical stability.²⁶ It possesses two different phases [Fig. 1(a)]. The room-temperature tetragonal phase (t-LLZO) with fully ordered Li-sublattice, leading to a low ion conductivity of ∼1 S cm⁻¹ K⁻¹ (ref. 27) in which ion migrations can only be promoted through phonon vibrations.²⁸ Above ∼900 K,^29,30 LLZO transforms into its cubic phase (c-LLZO) with a disordered Li-sublattice, in which the presence of Li vacancies significantly enhances the ion conductivity by four orders-of-magnitudes.²⁷ The cubic phase with Li-site disorder is the one that will be investigated in this work. Moving forward [Fig. 1(b)], the thermodynamic stability of c-LLZO can be further affected by the contents of Li⁺ (and oxygen stoichiometries associated with maintaining charge balance).^31,32

Fig. 1 Hierarchical sampling challenge in computational thermodynamics for disordered solids. Using garnet LLZO as an example: (a) the material possesses two (cubic and tetragonal) phases,^29,30 and in each phase (b) there is a variation in the chemical compositions of lithium. At any one composition, the Li⁺ ions can arrange themselves in a different combination of lattice sites, leading to a range of different atomic configurations of LLZO, each with its own enthalpy that is represented by the dots on the energy landscape. The line joining the lowest energy points across all the compositions represents the stability convex hull. More specifically, (c) different configurations of c-LLZO at a given Li content can be categorised based on the relative ratios of 24d to 96h site occupancies in the lattice. For a given ratio of Li site occupancies, (d) the PES of LLZO can be represented as multiple energy wells (grey curves), one for each configuration σ⁽ⁱ⁾. Typically, most of the current research focuses on computing the 0 K enthalpies, which correspond to the transformed PES³³ that is determined by the energies of the local minima H⁽ⁱ⁾ as represented by the solid horizontal lines. More accurate computation of materials' stability at finite temperature should also include (e) the vibrational contribution by replacing H⁽ⁱ⁾ with its corresponding temperature-dependent Helmholtz free-energy F⁽ⁱ⁾_vib(T).

Simply focusing on c-LLZO with full Li-stoichiometry, the disorder in the material originates from the total number (120) of lattice sites that can be occupied by the Li⁺ ions being much larger than the maximum number (56) of Li⁺ ions per unit cell of LLZO, resulting in an extremely large number³⁴ (∼7 × 10³⁴) of possible ways to arrange Li⁺ in the unit cell (commonly referred to as configurations). This means that configurational entropy could play a significant role in stabilising c-LLZO at high temperatures, and from a pure statistical thermodynamic perspective, the enthalpies for all these configurations need to be computed to obtain the converged configurational entropy³⁵ for c-LLZO with full Li-stoichiometry, which is computationally infeasible anyway. Crystallographically, different configurations of c-LLZO can be grouped based on the occupancy ratios of two types of Li-sites, namely the 24d and 96h‡ [Fig. 1(c)] (also see Fig. S1 in the ESI†). Importantly, Holland et al.³⁴ established a methodology to generate symmetrically unique configurations of c-LLZO with any site occupancy ratio, particularly in the range with optimum Li–Li distances to promote fast ion conduction. This not only provides the computational community with an exhaustive list of initial atomistic structures of c-LLZO, but also significantly reduces the configurational space, making it more attemptable to compute the thermodynamic stability of LLZO with the explicit inclusion of configurational entropy.

This is the first motivation behind the current work, whereby we examine the temperature-dependent thermodynamic stabilities for a subset of c-LLZO configurations with the 24d : 96h site occupancy ratios of 12 : 44 and 11 : 45, containing 4162 and 65 958 symmetry unique structures,³⁴ respectively. As in the standard structure prediction workflow, each configuration is first structurally optimised to its corresponding local energy minimum, which essentially provides us with the information of transformed PES at 0 K [Fig. 1(d)]. However, the thermodynamic quantity that truly dictates the stability of each configuration at finite temperature is its Helmholtz free energy,^7,20 which takes into account the contribution from lattice vibrations towards structural stability [Fig. 1(e)]. This can be computed with either a finite-displacement phonon³⁶ or molecular dynamics (MD) approach, but the high computational costs associated with these approaches, particularly when combined with a DFT energy model, prohibit these calculations from being individually performed for every configuration. Alternatively, the Helmholtz free energy of a representative configuration (typically the lowest energy one)²¹ may be used as a constant correction factor across the entire ensemble. The latter assumes Helmholtz free energies vary insignificantly amongst the configurations,²⁰ which may or may not hold.

The promising applications of LLZO as a solid electrolyte have already sparked significant interest in applying classical forcefields for conducting molecular dynamics investigations to understand its phase transition and Li-ion diffusion behaviours in pristine and aliovalent-doped LLZO, which have led to some important insights into the subtle interplays between its structural complexity and ion conductivities.^37–47 One of the key findings amongst these studies is the importance of increasing configurational entropy as reflected by Li redistributions in the LLZO lattice upon phase transition from the tetragonal to cubic phase.^41,42 Although both Ga and Al doping can stabilise LLZO in its high-conductive cubic phase, the Ga-substituted LLZO yields higher conductivity (10⁻³ S cm⁻¹) compared to the Al-substituted one (10⁻⁴ S cm⁻¹),⁴⁸ which may be attributed to the less repulsive nature of the Al³⁺ dopant that immobilises the nearby Li-vacancy and blocking the Li-diffusion pathway, as revealed from the forcefield simulation.^39,40 Technically, most of these classical forcefields are based on a Buckingham potential and point charge model. While this type of forcefield is fast to evaluate, Chen et al.⁴² pointed out the difficulty in pinpointing the exact Ta concentration and transition temperature to stabilise the cubic LLZO with Ta doping, due to the inherent limitation of the Buckingham potential. In this regard, Dai et al.⁴⁴ compared the performances of different forcefields for simulating fast ion conductors, and it is unsurprising to see the polarisable forcefield with an induced dipole model that captures the chemical-environment-dependent polarisation effects performs the best, but the development of an accurate polarisable forcefield is generally more challenging to non-specialists. Moreover, classical forcefields are typically fitted to the vicinity of local minimum, while this is beneficial for modelling the phononic properties, a good description of the long-ranged topological features of the PES is required for studying Li ion diffusions. Lastly, the pair potential also imposed limited flexibility to capture the important many-body atom–atom interactions in a complex system.⁴⁶

During the past decade, we have seen rapid progress in the development of machine learned forcefields (MLFFs)^49–64 with the emergence of the foundational model,⁶⁵ which offers accuracies in par with modern DFT but few orders of magnitude increase in the computational speed while overcoming the aforementioned limitations in classical forcefields, thus allowing us to bridge the gaps in length and time scales, as well as the compositional complexities in atomistic modelling. A few recent studies applied MLFF to further investigate the phase and ion diffusion behaviours in LLZO,^24,66–68 but critical quantitative insights into the role of configurational entropy in determining the thermodynamic properties of LLZO are still lacking. This makes it timely for us to harness the advantages of MLFF to tackle the sampling challenges and determine the temperature-dependent stabilities of LLZO, with the vibrational contributions to the configurational entropy accounted for explicitly, by following the exact statistical mechanics protocol.²⁰ This is the second motivation for the present study. More specifically, we will be using the SO3KRATES⁶² model to train our material-specific MLFF. It represents one particular implementation of the geometric graph neural networks⁶⁹ that use the message-passing scheme^59,70 to learn the information about connectivities in chemical structures. The architecture is made aware of the rotationally equivariant features by incorporating geometric priors based on the SO(3) rotational group -(special orthogonal group in three dimensions),^51,58,71,72 to improve the data efficiency, generalisability and stability in model training and atomistic simulations.^73,74 Furthermore, the network is combined with the self-attention mechanism⁷⁵ to bypass the need for full SO(3) tensor convolution in information exchange to increase the computational speed.

2 Theory

2.1 Thermodynamics of disordered solids

In this section, we briefly outline the fundamental thermodynamics of disordered solids. We begin with a set of structural models for the material of interest, each with symmetrically unique atomic arrangements. This set of structures are most commonly referred to as the configurations which will be collectively denoted as σ = {σ⁽¹⁾, σ⁽²⁾,…, σ⁽ⁿ⁾}. As the standard practice in (nearly all) atomistic modelling, we relax each configuration to its local energy minimum, giving us its corresponding enthalpy H⁽ⁱ⁾ at 0 K [Fig. 1(d)]. If we stop here, then the system is said to be stabilised by configurational entropy at T when the temperature-dependent configurational free-energy

G_conf(T) = H − TS_conf (1)

is negative, where the configurational entropy follows the Boltzmann relationship S_conf = k_B ln Z_conf. Here k_B is the Boltzmann constant, and the configurational partition functionwith β = 1/k_BT.§ Physically, Z_conf is a probabilistic measure of the likelihood that the system will simultaneously take the forms of all n distinct configurations at T. The system is more likely to remain disordered when Z_conf is large, meaning that all the configurations are energetically highly similar to each other, thus the system is not preferentially stabilised by a particular state with very low enthalpy. In this regard, what will be the enthalpy that goes into eqn (1). A sensible choice would be to take the temperature average for all configurationsrepresenting the enthalpy of the most likely configuration the system will take at a given T.

To further account for the vibrational contribution to the thermodynamic stability of each configuration under finite temperature [Fig. 1(e)],²⁰H⁽ⁱ⁾ in eqn (1) needs to be replaced by the corresponding vibrational (Helmholtz) free-energy F⁽ⁱ⁾_vib(T) = H⁽ⁱ⁾_T − TS⁽ⁱ⁾_vib, in which H⁽ⁱ⁾_T includes the internal energies at 0 K and those associated with thermal motion.⁷⁶ Here, we take H⁽ⁱ⁾_T = 〈U⁽ⁱ⁾(t)〉_MD, which is the time-averaging of the internal energy over an MD trajectory. With this, the configurational partition function becomes

and finally, the configurational free-energy

G_conf+vib(T) = F_vib(T) − TS_conf+vib, (5)

in which S_conf+vib = k_B ln Z_conf+vib. Physically, eqn (4) is equivalent of re-normalising the 0 K PES to its corresponding Helmholtz free-energy surface at a finite temperature.

Practically, there exist a multitude of approaches to compute F_vib.⁷⁶ The most common approach³⁶ is to first perform a frozen phonon calculation to determine the spectrum of phonon eigenfrequencies {ω⁽ⁱ⁾_λ} (where λ = (q,n) with q being the wavevector and n the band index), by diagonalising the Hessian matrix for the i-th configuration. Within the harmonic approximation, the vibrational free energy (ignoring H⁽ⁱ⁾_T) takes the following discrete sum:

Alternatively, one could take the integral form of the above expression and compute the vibrational free energy from the corresponding phonon density-of-states (PDOS) g⁽ⁱ⁾(ω) as^76–78

where N is the number of atoms in the simulation cell. In this work, we extract the PDOS from the Fourier transformation of the velocity autocorrelation function (VACF) that is computed from MLFF-MD simulation at the target temperature T aswith

There exist some notable disadvantages of computing F_vibvia VACF from MD simulations, which include: (a) MD should be performed at different temperatures T to extract the corresponding g(ω,T), if one is interested in the temperature-dependent Helmholtz free energies F_vib(T), for example, to determine the phase-transition temperature. (b) Large supercells, typically with ∼10³ atoms or more, need to be used and perform MD runs in sufficiently long times in principle, to obtain fully converged PDOS.⁷⁶ Both (a) and (b) translate to extremely high computational costs that are only achievable with forcefields, such as the MLFF adopted in this work. Finally, (c) the usage of eqn (7), which was established based on the harmonic approximation may not hold for strongly anharmonic systems. Nevertheless, it offers the advantage of accounting for the variations in PDOS that originate from the multiphonon interactions. Practically, running MD simulations with MLFF takes the advantage from the MD engines in the API (application programming interfaces) of many latest MLFF packages for atomistic simulations, that is either individually implemented or using other well-established APIs such as the atomic simulation environment (ASE),⁷⁹ as running long-time MD simulations is one of the key driving forces behind these MLFFs.

2.2 The SO3KRATES architecture

A neural network learns the contributions from individual atoms (E_i) towards the total energy of a structure with N atoms (as ) by exchanging the chemical and structural information amongst the constituent atoms through multiple interconnected hidden layers. Different neural networks vary in terms of how such information is exchanged while preserving the important symmetries in the chemical system. More specifically, this information exchange in SO3KRATES (see Fig. S2 ESI† for a detailed illustration) is achieved through T layers of Euclidean transformer blocks, which implements both the message-passing^69,70 and self-attention⁷⁵ mechanisms to iteratively update the per-atom feature asfrom which the total energy of the system is learned using a multi-layered perceptron (MLP) with the aggregated atomic features F^(T) = [f^(T)₁, f^(T)₂,…, f^(T)_N] as the input, i.e. E = MLP [F^(T)]. Here, the message m_ij holds the information about the j-th atom from all neighbours of i that is identified by a pre-defined cut-off radius r_cut. Before the first iteration, the feature vector is initialised purely with the chemical information (in this case, the atomic number z_i) of each atom through an embedding function as f⁽⁰⁾_i = f_emb(z_i) that maps z_i into a F-dimensional vector in the feature space, i.e.. As such, f encodes the invariant information of the system.

In addition, it is crucial for the SO3KRATES network to learn the geometric information of the chemical structure, that also contributes towards E_i and transforms equivariantly with respect to 3D rotations and translations. Here, the spherical harmonic coordinates (SPHC), the basis for the SO(3) rotational group, is used as an inductive basis to preserve the equivariance throughout the network. This is implemented, firstly by embedding the Cartesian vector r_ij = r_i − r_j for into its corresponding Euclidean variable (EV) as a multipolar expansion with

In the above equation, is the average number of atomic neighbours over the whole training data. ϕ_cut(||r_ij||) is a cosine cutoff function (eqn (17) in Frank et al.⁷⁴). r̂_ij = r_ij/||r_ij|| is the unit vector, and Y^m_l denotes the spherical harmonic function of degree with m ∈ [−l, −l + 1,…, l − 1, l]. Within each Euclidean transformer block, this SPHC is also updated per-degree-wise, similar to the feature vector as:

The messages [eqn (9) and (11)] are constructed using the multihead self-attention approach⁷⁵ (with the number of heads h being a network hyperparameter), to avoid the necessity of performing full SO(3) convolutions repetitively for L times along the feature vector dimension to exchange the invariant and equivariant information,⁷¹ bringing down the computational costs from to . More specifically, we write

m_ij = ϕ_cut(||r_ij||)α_ijf_j, (12)

for the feature message and

m_ij,lm = ϕ_cut(||r_ij||)α^l_ijY_l^m(r̂_ij), (13)

for the EV. Here, the self-attention coefficient is calculated as

α_ij = k_i^T(w_ij⊙q_j), (14)

in which the query and key vectors are computed as q_j = Qf_j and k_i = Kf_i with Q and K being learnable matrices of type . The weight matrix w is also learnable as the sum of the outcomes from two MLPs as

in which contracts each degree in c_ij to the zeroth degree (eqn (14) in Frank et al.⁷⁴). This completes the Euclidean attention block in the Euclidean transformer. The outputs from the attention block (f^attn_i and χ^attn_i) are subsequently passed into the interaction block, which provides the network with additional degree-of-freedoms to exchange the invariant and equivariant information. The outcomes are the refinement df_i and dχ_i, which are to be combined with the initial inputs to generate the final outputs for a given Euclidean transformer block:

f^(t+1)_i = f^(t)_i + df_i,

and

χ^(t+1)_i = χ^(t)_i + dχ_i.

Finally, autodifferentiation is employed to calculate the atomic forces to ensure energy conservation.⁵⁶

3 Results

3.1 Accuracy of the model for structural optimisation

Using the equivariant features as an inductive basis in the neural network has demonstrated its advantages in achieving outstanding model performance and data efficiency.^62,71 When training the SO3KRATES model on the structural optimisation trajectories using the 12_44_opt set (see Section 5.2 for details), it is found from the learning curves (Fig. S3, ESI†) that using as little as 2000 data points (corresponding to just 2% of the database) is sufficient for the model to reach the meV-level of accuracy in predicting the energies and forces. Beyond that, the model performances only improve slightly with further increase in the size of the training set, and in the subsequent discussions, the model that is applied for structural optimisations is trained with 4000 data points. The accuracy of this particular MLFF is further demonstrated in Fig. 2 which shows both the energies and atomic forces predicted by the MLFF correlate very well with those calculated from DFT, reaching an error level of 7.48 meV and 4.87 meV per (Å atom) in MAE (mean-absolute-error) with 9.62 meV and 6.72 meV per (Å atom) in RMSE (root-mean-squared-error) for energy and force predictions, respectively.

Fig. 2 Accuracy of the MLFF trained from structure optimisation trajectories. Comparison of the MLFF-predicted and the target DFT (left) energies and (right) atomic forces across 700 randomly selected configurations from the 12_44_opt set. The model is trained with 4000 data.

To gain some insights into how the performance of the MLFF is affected by the model flexibility, we trained extra versions of MLFFs by changing the dimensionality of the EV χ_i through varying the degrees of spherical harmonics l. The resulting model performances, which are compared in Table 1, show that compared to the atomic forces, the accuracy in predicting energy is more sensitive to the degrees of multipolar expansions used for constructing the EVs. This can be due to the fact that the models are trained using a loss function that is strongly biased towards accurate force predictions. The largest improvement in the model accuracy can be seen when the dipole moment (l = 1) is incorporated in additional to the monopole moment (l = 0), signifying the importance of learning the orientational information in the chemical environments. Nevertheless, the relative contributions from higher order multipoles reduce as they are progressively introduced into the model. Interestingly, EVs that are expanded to the hexadecapoles while ignoring the monopole contribution shows the best accuracy in energy prediction, with only a small deterioration in forces (compared to the case where the monopole is included). We speculate that this may be because the structural framework of LLZO is made of multiple LiO₄, LiO₆, LaO₆ and ZrO₆ polyhedra, which involve many chemical bondings that are strongly directional. In this case, including the monopole term, which is isotropic in nature, is not beneficial to capture the directional bonding characteristics, causing the model to become overfit resulting in worse performance. As such, this choice of l is used throughout the rest of this study.

Table 1 Effect of varying the degrees of spherical harmonics l applied for constructing the EVs on the performances of the MLFF models, which are tested on 700 randomly selected data from the 12_44_opt set

l	Energy [meV]		Force [meV per (Å atom)]
	MAE	RMSE	MAE	RMSE
a This corresponds to the model shown in Fig. 2.
{0}	25.0	31.1	8.68	12.3
{0,1}	15.0	19.0	6.36	8.70
{0,1,2}	11.7	14.9	4.92	6.72
{0,1,2,3}	10.3	14.2	4.45	6.14
{1,2,3}^a	7.48	9.62	4.87	6.72

---

3.2 Accuracy of the model for molecular dynamics simulations

Unfortunately, the MLFF trained on the structural optimisation trajectories is not sufficiently extrapolative to produce stable MD runs,⁷³ because the latter can sample a much wider configurational space with broader distributions in energies and atomic forces that are far beyond the training data covered in the optimisation trajectories. Hence a separate MLFF is trained on a single AIMD (ab initio molecular dynamics) trajectory which heats up a c-LLZO structure from 10 to 2100 K.⁶³ The lattice constants of c-LLZO were kept fixed during the simulation as the examination of phase transition behaviours in LLZO is not the main goal of the current study.²⁴ Similar to the training of MLFF for structural optimisations discussed in Section 3.1, the top row of Fig. 3 again shows the model is capable of reaching an meV-level of accuracy with a small number of 1000 training data.

Fig. 3 Accuracy of the MLFF trained for molecular dynamics simulations. Correlations for energies (a, c) and atomic forces (b, d) between those predicted from the MLFF and the respective target values from DFT. The top row (a, b) demonstrates the model accuracy using all 5200 data from the AIMD trajectory used for training the model, whereas the bottom row (c, d) demonstrates the model transferability with data from 50 AIMD trajectories under the NVT ensemble that have not been used for model training.

Apart from the model being accurate, the transferability of the MLFF is also critical, if we would like to use this particular MLFF to determine the vibrational contributions to the configurational entropies from all other configurations in the set. For this, at the bottom row of Fig. 3, we show that the model is capable of reproducing the DFT energies and forces from the 300 K AIMD trajectories of another 50 configurations that are randomly selected from the 12_44_opt set. Note that the prediction errors for this set are seemingly lower than those computed on the AIMD trajectory for model training, because the latter covers much wider ranges in energies and forces. In addition to the equivariant architecture of the SO3KRATES network, the excellent transferability of the trained MLFF could also be attributed to the fact that, despite the configurational space for c-LLZO being very large, the short-ranged chemical environments,⁵¹ which are expected to play a dominant contribution to the cohesive energy of a crystal, are much less diverse across configurations.

We also tested the stabilities and performances of using the as-trained MLFF in molecular dynamics simulations across the temperature range (300–1500 K) that we are interested in. Because the key quantity that is needed for computing F_vib is the PDOS [eqn (7)], in Fig. 4, we compare the temperature-dependent PDOS for one configuration of c-LLZO that is computed from AIMD and MLFF-MD via the VACF, respectively. From Fig. 4, it can be seen that the overall shapes of the PDOS from AIMD are well reproduced by the MLFF-MD, particularly the highest PDOS peak that is centred at ∼40 meV. In the phonon energy range of 50–70 meV, MLFF tends to give slightly lower (higher) PDOS values at the lower (higher) temperatures, compared to the ones computed from AIMD. Recent investigation in Li₃YCl₆⁸⁰ shows the appearance of low-frequency boson-like peaks⁸¹ in its PDOS as a signature of strong vibrational anharmonicity and coupled vibration between Li and the rest of the crystal framework plays a key role towards the emergence of the superionic state in this material. However, we did observe similar behaviour from the PDOS of LLZO. Finally, no significant shift in the PDOS peak(s) is observed in both PDOS spectra across the temperature range investigated, suggesting that thermally induced phonon renormalisation in c-LLZO is not expected to be significant. Nevertheless, more investigations of the phononic properties of c-LLZO are subject to future studies.^82,83 In Table 2, we further compare the key vibrational thermodynamic quantities that are extracted from AIMD and MLFF-MD simulations. Statistically, compared with DFT over the entire temperature range investigated, the ML potential is found to be slightly more attractive by 121.5 meV per atom (as measured by the time-averaged internal energy 〈U(t)〉_t), which may be attributed to the higher PDOS at the higher energy range produced from MLFF-MD as shown in Fig. 4. At the same time, the vibrational entropy S_vib computed from MLFF-MD is slightly larger by 18.6 meV per (atom K) compared to the DFT value. It is also shown in Table 2 that, at around 1500 K where c-LLZO is the more stable phase, the vibrational contribution (TS_vib) to the free energy, albeit smaller, is on the same order of magnitude as the internal energy 〈U(t)〉_t, signifying the importance of accounting for the vibrational contribution in the accurate determination of the material stability.

Fig. 4 Performance of the MLFF in computing the vibrational properties of LLZO. Comparison of phonon DOS computed for one configuration in the 12_44_opt set at different temperatures from (a) AIMD and (b) MLFF-MD.

Table 2 Comparison of the internal energies (〈U(t)〉_t), vibrational entropies (S_vib) and free energies (TS_vib) computed from the MD trajectories at different temperatures using DFT and the trained MLFF for one c-LLZO configuration

Temp. [K]	〈U(t)〉^DFTt [eV per atom]	〈U(t)〉^MLFFt [eV per atom]	S ^DFTvib [meV per (atom K)]	S ^MLFFvib [meV per (atom K)]	TS ^DFTvib [meV per atom]	TS ^MLFFvib [meV per atom]
300	−7.802	−7.820	0.912	0.928	273.65	278.46
500	−7.748	−7.808	0.905	0.935	452.58	476.61
700	−7.709	−7.795	0.893	0.934	625.45	653.90
900	−7.672	−7.782	0.896	0.918	806.52	852.91
1000	−7.656	−7.776	0.902	0.912	902.36	911.69
1100	−7.636	−7.769	0.887	0.910	975.64	1001.3
1200	−7.625	−7.763	0.901	0.904	1081.6	1084.9
1300	−7.594	−7.756	0.884	0.902	1149.1	1173.1
1400	−7.567	−7.750	0.886	0.897	1240.3	1256.1
1500	−7.538	−7.743	0.876	0.888	1312.4	1331.3

---

With the demonstrated accuracy of the trained MLFF, we are able to further examine the effect of supercell size on the PDOS. In Fig. S4 of the ESI,† we show that the changes in the computed PDOS using the supercell with up to 3840 atoms are again minimal, particularly with no sign of soft phonon modes emerging from the PDOS simulated with a larger supercell, thus suggesting that using a single unit cell of c-LLZO is sufficient to converge the PDOS calculations via VACF.

Although this is not the primary focus of the current study, in Section S2 of the ESI† we also compared the Li ion diffusivity and the corresponding activation energies that are further compared to those values extracted from previous experimental and theoretical results, which also showed good agreement between our MLFF-MD results and previously published data.

Finally, we highlight the computational efficiency delivered by MLFF. With AIMD, simulating a 4-ps-long trajectory using a single unit cell of c-LLZO on 56 CPUs over two Broadwell nodes takes ∼34 hours walltime. In contrast, using MLFF executed on a single V100 GPU, we are able to sample a 10-ps-long trajectory in ∼7 minutes, making the comprehensive thermodynamic sampling of disorder solids, such as LLZO, within the reachable computational capabilities.

3.3 Convergence of configurational free energy without the vibrational contribution

Before the trained MLFF can be applied to compute the configurational free energy of c-LLZO, its capability to reproduce the DFT local minima [cf.Fig. 1(d)] from structural optimisations needs to be further verified. Since our primary interest is the configurational free energy, instead of comparing the structural similarities between the DFT- and MLFF-optimised crystal structures of every c-LLZO configuration, we focus on the differences in energies of the local minima that are sampled by these two different energy models. More specifically, we compare the similarity in the distributions of the energies p(H), as the partition function [eqn (2)] can be re-written in the following integral form for a continuous energy spectrum:
and the result is shown in Fig. 5. Here, all the starting structures from the 12_44_opt have been re-optimised with MLFF to generate a new set of local minima on the MLFF-PES.

Fig. 5 Performance of MLFF for structure optimisations. (Solid lines) Distributions of energies p(E) for all the local minima on the PES for the 12_44_opt set. These distributions are shown as the Gaussian Kernel Density Estimates (KDEs) of the underlying data, and KDE is used to present other similar data throughout the rest of this paper. (Dashed lines) Evolutions of the configurational entropic contribution to the free energy (−TS_conf) at 300 K as more configurations are included (in the order of increasing energy) in the summation [cf.eqn (2)].

It can be seen from Fig. 5 that, firstly, the overall shapes of p^MLFF(H) and p^DFT(H) match well with each other, which again demonstrates the trained MLFF for structural optimisation is very accurate. Because structural optimisation essentially traces the shape of the PES near every local minimum, it means this topological feature has been well reproduced by the MLFF, even though it is trained with only 4% of the available data. Secondly, the energies of the local minima tend to be slightly underestimated by MLFF. This could be due to a different structural optimiser and looser convergence criteria that were used in the structural optimisation with MLFF. Nevertheless, as shown by the dashed lines in Fig. 5, the configurational entropies computed from DFT and MLFF converge to the same value as soon as the highest peak in p(E) is passed, and Table 3 further shows that all thermodynamic quantities computed from MLFF are essentially identical to those from DFT.

Table 3 Comparison of the configurational entropies and free energies computed for all c-LLZO configurations in the 12_44_ opt set using DFT and MLFF at 300 and 1500 K

Temp. [K]		〈H〉T [eV per atom]	S conf [meV per (atom K)]	−TSconf [eV per atom]	G conf [eV per atom]
300	DFT	−7.839	0.7112	−0.2134	−8.052
	MLFF	−7.839	0.7118	−0.2135	−8.053
1500	DFT	−7.839	0.7167	−1.075	−8.914
	MLFF	−7.839	0.7168	−1.075	−8.914

---

For comparison, we trained a separate MLFF for t-LLZO and computed its vibrational free energies (see Table S3 of the ESI,† and because this structure is not disordered, no configurational contribution needs to be included). It can be seen from Table 3 that, at 1500 K, although disorder has contributed an additional 1 eV per atom to the stability of c-LLZO leading to a configurational free energy of −8.914 eV per atom, it is still larger than the vibrational free energy of −9.147 eV per atom for t-LLZO.

3.4 Vibrational contribution to the configurational entropy

Now, we are in the position to examine how thermal vibrations renormalise the PES into the corresponding temperature-dependent free energy surfaces, and subsequently affect the configurational contribution to the thermodynamic stability of site-disordered c-LLZO. We begin with a brief look on how PDOS varies across the configurational space, and in Fig. 6 we have plotted the temperature-dependent configurational-averaged PDOS⁸⁴ for the 12_44_opt set, which are computed as
with the corresponding standard deviation
in which with p̃⁽ⁱ⁾(T) = exp(−〈U⁽ⁱ⁾(t)〉_t/k_BT). It can be seen from Fig. 6 that, at a given temperature, the variation in the PDOS as measured by σ[g(ω)] is substantial and does not change significantly across the entire phonon energy range, meaning the vibrational behaviour of the LLZO crystal lattice is very sensitive to how the Li ions are arranged in the crystal. This supports the necessity to compute the vibrational free energy per configuration individually as opposed to approximating it with the value of a single configuration for very accurate results. Nevertheless, σ[g(ω)] are larger at 300 K compared to those at 1500 K, possibly because the latter is the temperature at which LLZO is stable at the cubic phase. Fig. 6 also clearly indicates the increase of PDOS in the energy range above 50 meV from 300 to 1500 K, reflecting the excitation of high energy phonons at the elevated temperatures.

Fig. 6 Configurational dependence of phonon behaviours PDOS that is computed as the weighted average of PDOS for all c-LLZO configurations with an Li site occupancy of 12 : 44 [Section 3.4 and the associated standard deviations across the range of phonon energies covered]. The data are presented for calculations performed at (a) 300 K and (b) 1500 K, respectively.

With these, Fig. 7 further demonstrates the distributions of vibrational free energies computed for all configurations in the 12_44_opt set. Compared with the distributions of 0 K-enthalpies Fig. 5, it shows that as the temperature increases, the (free) energy distributions gradually widened and shifted to lower energy ranges, indicating the structure is progressively stabilised by lattice vibrations. The changes in the shapes of the (free) energy distributions further show that the vibrational entropic contribution to the free energy is not a simple constant correction factor across the configurational space. By integrating the free energy distributions, the configurational entropy and free energy that include the vibrational contribution at 300 and 1500 K can be obtained and are listed in Table 4. Using the free energies (G_conf+vib) determined in this way, the cubic phase of LLZO is now correctly predicted to be more stable than the tetragonal one by 973 meV at 1500 K (as opposed to being 233 meV less stable with G_conf only), confirming the importance of lattice vibrations in stabilising the cubic phase at high temperature. In contrary, at 300 K, the cubic phase is also predicted to be more stable than the tetragonal phase by 106 meV based on G_conf+vib, and the later was predicted to be more stable by 147 meV if the vibrational contributions were omitted for c-LLZO. This could be indicating that the lattice vibrations and configurational disorders are strongly competing in determining the phase stability at 300 K, which is worth further investigating in the future.

Fig. 7 Vibrational free energy landscapes Probability distributions p(F_vib) (solid lines) of the vibrational free energies F_vib for all c-LLZO configurations with an Li site occupancy of 12 : 44 at 1500 and 300 K. Note that the two figures are arranged from left to right in the order of increasing F_vib. The dashed line shows the convergence of the entropic (vibrational and configurational) contribution to the free energy.

Table 4 Configurational entropies and free energies, with vibrational contributions taken into account, for all c-LLZO configurations with an Li site occupancy of 12 : 44. The expectation value of the configurational averaged free energy at a given temperature, 〈F_vib〉_T, is calculated using eqn(3) by replacing the enthalpy H⁽ⁱ⁾ with the vibrational free energy F⁽ⁱ⁾_vib [eqn(7)] for each configuration, such as, G_conf+vib = 〈F_vib〉_T − TS_conf+vib

Temp. [K]	〈Fvib〉T [eV per atom]	S conf+vib [meV per (atom K)]	−TSconf+vib [eV per atom]	G conf+vib [eV per atom]
300	−8.107	0.6578	−0.1973	−8.305
1500	−9.093	0.6863	−1.029	−10.12

---

3.5 Model transferability to c-LLZO with a different Li site occupancy

Having demonstrated the robustness of the MLFFs that are trained on the structural optimisation and AIMD trajectories for c-LLZO structures with an Li site occupancy of 12 : 44, in this section, we shall apply these MLFFs to explore the PES of a larger set of c-LLZO with the Li site occupancy of 11 : 45 that contains a total of 65 958 unique configurations.³⁴ Similar to Fig. 5, in Fig. 8(a), we first compare the distribution of the local minimum energies p(E) that are obtained through structural optimisation using MLFF, using a subset of c-LLZO containing 19 333 configurations, for which the corresponding DFT results are available. Once again, we observe good agreement between p^MLFF(E) and p^DFT(E), whereby the key features of the peak and shoulders in p^DFT(E) are well reproduced in p^MLFF(E). This suggests that the as-trained MLFF can be well transferred to model c-LLZO of a different Li site occupancy, and again demonstrates the dominant contributions to the cohesive energies of the crystal is the short-range chemical environment rather than the overall site-occupancy in the unit cell.

Fig. 8 Energy landscape for cubic LLZO with an Li site occupancy of 11 : 45. (a) Benchmarking the performance of the MLFF trained on the 12_44_opt for optimising a subset of c-LLZO structures (19 333 data points) with an Li site occupancy of 11 : 45. The solid lines compare the distribution of energies p(E) for the local minima on the PES as obtained from DFT and MLFF, whereas the dashed lines show the convergence of configurational entropic contribution to the free energy at 300 K. (b)–(d) Distributions of enthalpies (b) and vibrational free energies (c) and (d) for all 65 958 configurations of c-LLZO with an Li site occupancy of 11 : 45 at 0, 300 and 1500 K, as computed using MLFF, respectively, and the convergence of the configurational entropic contribution to the free energies calculated from the corresponding energy distributions. In (b), the DFT data [as shown in (a)] are also included for comparison.

The most significant advantage offered by MLFF is its speed that enables us to perform the structural optimisation, together with MD samplings at 300 and 1500 K for the entire set of 65 958 c-LLZO configurations in ∼344 days on a single V100 GPU. This is much less than the ∼810 days that it cost us just to perform structural optimisations on 19 333 configurations (29% of the dataset) using periodic DFT. The results are presented in Fig. 8(b–d) and Table 5.

Table 5 Configurational entropies and free energies for all c-LLZO configurations with an Li site occupancy of 11 : 45

Without vibrational contribution
Temp. [K]	〈H〉T [eV per atom]	S conf [meV per (atom K)]	−TSconf [eV per atom]	G conf [eV per atom]
300	−7.839	0.9494	−0.2848	−8.123
1500	−7.839	0.9549	−1.432	−9.270

---

With vibrational contribution
Temp. [K]	〈Fvib〉T [eV per atom]	S conf+vib [meV per (atom K)]	−TSconf+vib [eV per atom]	G conf+vib [eV per atom]
300	−8.106	0.8792	−0.2638	−8.369
1500	−9.092	0.9195	−1.379	−10.47

---

The p(H) curves in Fig. 8(b) show that the full set [as given by p^MLFF(H)] contains more lower energy configurations than the DFT set [as given by p^DFT(H)]. Based on these two distributions, at 300 K, the computed configurational entropic contributions to the free energy −TS_config is −0.2528 eV per atom from the DFT set, which is 11.2% lower than the converged value of −0.2848 eV per atom from MLFF sampling [also see Table 5]. This level of error is qualitatively consistent when ∼20k of MLFF-optimised local minimum energies are used to determine −TS_config [see Fig. S8, ESI†]. This underestimation in the configurational entropy can be primarily due to incomplete sampling, and this can be tested by using the KDE (Kernel Density Estimate) density of p^DFT (H) to sample a full set of 65 958 enthalpy values and compute the corresponding configurational entropy at 300 K, which gives an estimated −TS_config of −0.2843 eV per atom or 0.176% relative error. Note that the size of the DFT-optimised set (19 333) was not chosen based on any statistical insight but merely an arbitrary point at which we thought there would be sufficient data for benchmarking purposes and decided to stop carrying on with more DFT calculations. In the above analysis, as shown by the density estimations, the DFT set is definitely oversampled, consuming more computing time than it is necessarily needed. A promising future direction is to explore generative models to enhance the sampling efficiency for disordered solids,^85–87 but it is beyond the scope of the present work.

Moving on to the vibrational free energy landscapes at finite temperatures, Fig. 8(c and d) show again that the distributions of F_vib widen and shift toward the lower energy range as the temperature is increased from 300 to 1500 K, which is consistent with the observations from Fig. 7. On the other hand, although the number of configurations in this set of c-LLZO is an order of magnitude more than the other set with an Li site occupancy of 12 : 44, the converged G_conf+vib values for this set of c-LLZO are only 64 and 350 meV lower than those for the 12 : 44 set at 300 and 1500 K, respectively (see Table 5 and Table 4). This is because the distributions of F_vib [Fig. 8(c and d)] are much wider than the thermal energy k_BT at the respective temperatures, such that many high-energy metastable configurations contributed insignificantly towards the partition function.

A natural extension from here on will be to apply the as-trained MLFFs to compute the configurational entropies for c-LLZO with other Li site occupancies [Fig. 1(c)], from which the total configurational entropy across the global PES may be obtained from the superposition of energy basins.^88–91 However, this is merely an exercise of exhausting computational powers, and more informative sampling strategies should be explored to achieve this goal efficiently to advance the field of computational thermodynamics of disordered solids.

4. Conclusions and future outlooks

In this work, for the first time, we are able to follow the exact statistic thermodynamic protocol, to computationally determine the combined contributions from vibrational and configurational entropies towards the phase stabilities of cubic garnet electrolyte LLZO with two specific Li site occupancies. This is enabled by using the equivariant neural network SO3KRATES that is built upon the transformer architecture, delivering high data efficiency, regression accuracy, enhanced stability in MD simulations, as well as good model transferability. Our results show uniquely, that the configurational free energies determined from the 0 K-enthalpies, which correspond to the energies of the local minima on the PES, are insufficient to correctly order the cubic LLZO over its tetragonal counterpart as the stable phase at 1500 K. This can be remedied instead, by computing the configurational free energies from the Helmholtz free energies, to account for the effect of lattice vibrations in renormalising the 0 K PES to the corresponding finite-temperature free energy surface. Most importantly, our results clearly demonstrate significant variations in the PDOSs, which translate to a non-uniform shift from H to F_vib, across different configurations, rendering the necessity to compute F_vib individually for every LLZO configuration. This task is certainly only feasible by using a fast and stable forcefield, like the one being developed in this work.

More specifically, we have chosen to compute the Helmholtz free energies by convoluting the PDOS from MLFF-MD simulations. While this offers the advantage to account for the anharmonicity in lattice vibrations, we cannot fully rule out the possibility of Li ions diffusing into other neighbouring sites due to their highly mobile nature. This means that F_vib calculated in this way may have additionally included contributions from neighbouring local minima that may belong to other Li site occupancies, especially at high temperatures. This could be checked, in the future, by recalculating F_vib based on quasiharmonic approximations from frozen-phonon calculations. Another avenue worth exploring is to train the MLFFs to directly predict F_vib using the data generated in this work, which can immediately offer an order of magnitude speed up in the sampling efficiency.

Despite our explicit account for the vibrational contributions to the phase stabilities of c-LLZO at finite temperatures, we have only considered LLZO with perfect stoichiometry and the effects such as dopants together with Li and O non-stoichiometries in lowering the stable temperature of c-LLZO have been neglected.^31,32,92,93 These factors not only diversify the chemical space⁹⁴ but also increase the number of configurations factorially, making the brute-force approach adopted in this work no longer feasible even with an efficient MLFF. The use of foundational models⁶⁵ trained across the chemical space of possible dopants in LLZO, that is further combined with generative models,^85–87 may offer the much needed solutions in the near future.

Finally, we would like to reiterate a key difference between the traditional atom–atom forcefield⁹⁵ and MLFF. In the language of artificial intelligence, the former can be considered as ‘rule-based’ models, where there is a clear physical motivation behind the design of each functional employed in the forcefield for computing the strengths of atom–atom interactions, such as the two-body bond stretching and dispersion interactions. In contrast, MLFF developed based on neural network architecture is largely relying on information exchange to surrogate a high-dimensional function, or an ansatz, that best describes the behaviours of the training data. In this regard, the development and adaptation of explainable artificial intelligence in computational chemistry and materials science will be an important avenue to explore for advancing our understanding of the structure–property relationships for functional materials at the atomistic level.⁹⁶

5 Methodologies

5.1 Ab initio calculations

All DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP)⁹⁷ with the standard projector augmented wave (PAW) method⁹⁸ and the PBE -(Perdew–Burke–Ernzerhof)⁹⁹ exchange–correlation functional. We employed the GW-version of the pseudopotential for each element at an energy cut-off of 500 eV. Geometry optimisations were performed only to relax the atomic positions while keeping the lattice constants unchanged, such that the latter are consistent across all training/validation data. The optimisations were terminated when the energy difference between two successive steps is below 10⁻⁴ eV. All calculations were performed without spin polarisation and only a single Γ-centered K-point was adopted in order to speed up the computation of a large set of structures. The K-point setting was chosen by considering that the LLZO structures are of a large unit cell with lattice constants a ∼12–14 Å.

To generate the data for training the MLFF for MD simulations, one structure of cubic LLZO with an Li-site occupancy of 12 : 44 was first randomly selected. Then, an AIMD simulation was carried out from 10 to 2610 K at 1 fs per step for 5.2 ps using the velocity-scaling algorithm to generate a total of 5200 frames for model training and validations. A single unit cell of LLZO is used for AIMD and its lattice constants were kept fixed during the simulation. To achieve a stable MD run using MLFF at a target temperature T_c, it is critical that the MLFF model can learn the topology of the PES over an extensive energy range, such that the model ‘knows’ how to restore to a stable structure if some overly distorted structures are generated during the MD run. This can be achieved by including structures that have been sampled over a wide temperature range in AIMD to generate the data to train the MLFF for MD simulations.⁷³ For benchmarking the phonon DOS using VACF, AIMD simulations were performed under the NVE ensemble using the velocity Verlet algorithm at the target temperatures at 1 fs per step for 4 ps. To further test the model transferability for running MD simulations on other structures in the 1244 opt set (see Section 5.2 for details), 50 structures were randomly selected and AIMD simulations were performed under the NVT ensemble at 320 K for 4 ps at 1 fs per step, using an Anderson thermostat with a collision probability of 0.5. Note that the AIMD trajectories generated under this ensemble are not suitable for computing the phonon DOS, because the atomic forces are ‘contaminated’ by the coupling to the heat bath. As such, they will only be used for comparing the accuracies of energies and atomic forces predicted by the MLFF model.

5.2 Training of the SO3KRATES networks

The SO3KRATES models are trained with the following settings. The dimension of the feature vector F is 132, and h = 4 heads are used in the message-passing updates. The cut-off distance r_cut = 5 Å is chosen for identifying neighbouring atoms. l = {1,2,3} is selected for constructing the EVs.

For the task of structural optimisations, the SO3KRATES model is trained with T = 2 layers. The 12_44_opt set is used for model training and validation. Briefly, all cubic LLZO structures with a 24d : 96h site occupancy ratio of 12 : 44 are retrieved from the dataset published by Holland et al.³⁴ Then the atomic positions in each structure are randomly disturbed using a method implemented in the pymatgen Python API,¹⁰⁰ before it is optimised to its nearby local minimum using DFT. The random distributions help to generate a more diversified set of atomic environments for fitting a robust MLFF model. In total, the 12_44_opt contains 95 271 data points, and 1000 are selected for model validation. Model training is terminated after 500 epochs.

For the task of molecular dynamics simulations, the SO3KRATES model is trained with T = 3 layers, with 1000 data points used for training and 1200 for validation. Model training is terminated after 1000 epochs.

All data points for training and validation were randomly selected from the full (structural optimisation or AIMD trajectories) set without applying any filtering procedure. We note that more informed and robust strategies to select training data exist, for example, in our previous work,¹⁰¹ structural descriptor and kernel-based ML combined with furthest-point sampling was applied to select a set of structurally most distinctive molecular crystals to enhance the model accuracy while reducing the amount of the training data used. However, using this approach would require the hyperparameters for the kernel to be separately optimised to achieve the best expressivity. Computationally, large amounts of disc storage and memories are also necessary to store the pre-computed structural descriptors and/or kernels, which can be particularly significant for LLZO with more than 100 atoms in the unit cell. This means separate and devoted efforts are required to implement this workflow. As such, we did not implement the furthest-point sampling strategy to optimise the training and validation sets in this work, but it is worth exploring in the future.

All models are trained by minimising the following loss function that combines the energy and atomic forces:

in which the trade-off parameter β = 0.99. The ADAM (adaptive moment estimation)¹⁰² optimiser is used with an initial learning rate of 10⁻³, which is decreased by a factor of 0.9 for every 10⁴ training steps. All model trainings are carried out with a single NVIDIA V100 GPU.

5.3 Structural optimisations and MD simulations with MLFF

Structural optimisations are performed using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm¹⁰³ to optimise the atomic positions, and are terminated when the maximum atomic force is smaller than 0.08 eV Å⁻¹. MD simulations at the target temperature are carried out under the NVE ensemble for 6 ps at 1fs/step using the velocity Verlet algorithm.

Data availability

The raw data for this work can be freely accessed from https://dataverse.harvard.edu/dataverse/unsw_sse_database. Python scripts used for data generations and analysis in this work can be found at https://github.com/yangjackie/futuremat_public.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

All authors would like to acknowledge the financial support from the Australian Research Council under the projects DP220103229, as well as the UNSW Materials and Manufacturing Futures Institute. Computations in this work are conducted on the Gadi supercomputer at NCI Australia that is supported by the UNSW Merit Allocation Scheme and the National Computational Merit Allocation Scheme (NCMAS) under the project dy3, as well as the local katana cluster supported by the Research Technology Service, UNSW.

References

1. L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid and M. G. Kanatzidis, Ultralow thermal conductivity and high thermoelectric figure of merit in SnSe crystals, Nature, 2014, 508, 373.
2. T. Zhu and E. Ertekin, Mixed phononic and non-phononic transport in hybrid lead halide perovskites: glass-crystal duality, dynamical disorder, and anharmonicity, Energy Environ. Sci., 2019, 12, 216.
3. M. M. Ugeda, A. J. Bradley, S.-F. Shi, F. H. Da Jornada, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, F. Wang, S. G. Louie and M. F. Crommie, Giant bandgap renormalization and excitonic effects in a monolayer transition metal dichalcogenide semiconductor, Nat. Mater., 2014, 13, 1091.
4. X. Jiang, Q. Zheng, Z. Lan, W. A. Saidi, X. Ren and J. Zhao, Real-time GW-BSE investigations on spin-valley exciton dynamics in monolayer transition metal dichalcogenide, Sci. Adv., 2021, 7, eabf3759.
5. P. You, D. Chen, X. Liu, C. Zhang, A. Selloni and S. Meng, Correlated electron-nuclear dynamics of photoinduced water dissociation on rutile TiO₂, Nat. Mater., 2024, 23, 1100.
6. A. Otero-De-La-Roza and E. R. Johnson, A benchmark for non-covalent interactions in solids, J. Chem. Phys., 2012, 137, 054103.
7. J. Nyman, O. S. Pundyke and G. M. Day, Accurate force fields and methods for modelling organic molecular crystals at finite temperatures, Phys. Chem. Chem. Phys., 2016, 18, 15828.
8. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, A look at the density functional theory zoo with the advanced GMTKN55 database for general main group thermochemistry, kinetics and noncovalent interactions, Phys. Chem. Chem. Phys., 2017, 19, 32184.
9. S. A. Tawfik, T. Gould, C. Stampfl and M. J. Ford, Evaluation of van der Waals density functionals for layered materials, Phys. Rev. Mater., 2018, 2, 034005.
10. Y. Zhang, D. A. Kitchaev, J. Yang, T. Chen, S. T. Dacek, R. A. Sarmiento-Pérez, M. A. Marques, H. Peng, G. Ceder, J. P. Perdew and J. Sun, Efficient first-principles prediction of solid stability: Towards chemical accuracy, npj Comput. Mater., 2018, 4, 9.
11. S. L. Price, Predicting crystal structures of organic compounds, Chem. Soc. Rev., 2014, 43, 2098.
12. D. H. Bowskill, I. J. Sugden, S. Konstantinopoulos, C. S. Adjiman and C. C. Pantelides, Crystal structure prediction methods for organic molecules: State of the art, Annu. Rev. Chem. Biomol. Eng., 2021, 12, 593.
13. A. R. Oganov, A. O. Lyakhov and M. Valle, How Evolutionary Crystal Structure Prediction Work, and Why, Acc. Chem. Res., 2011, 44, 227.
14. L. Wei, N. Fu, E. M. Siriwardane, W. Yang, S. S. Omee, R. Dong, R. Xin and J. Hu, TCSP: a template-based crystal structure prediction algorithm for materials discovery, Inorg. Chem., 2022, 61, 8431.
15. A. Pulido, L. Chen, T. Kaczorowski, D. Holden, M. A. Little, S. Y. Chong, B. J. Slater, D. P. McMahon, B. Bonillo, C. J. Stackhouse, A. Stephenson, C. M. Kane, R. Clowes, T. Hasell, A. Cooper and G. M. Day, Functional materials discovery using energy-structure-function maps, Nature, 2017, 543, 657.
16. C. Collins, M. Dyer, M. Pitcher, G. Whitehead, M. Zanella, P. Mandal, J. Claridge, G. Darling and M. Rosseinsky, Accelerated discovery of two crystal structure types in a complex inorganic phase field, Nature, 2017, 546, 280.
17. W. Sun, C. J. Bartel, E. Arca, S. R. Bauers, B. Matthews, B. Orvañanos, B.-R. Chen, M. F. Toney, L. T. Schelhas, W. Tumas, J. Tate, A. Zakutayev, S. Lany, A. M. Holder and G. Ceder, A map of the inorganic ternary metal nitrides, Nat. Mater., 2019, 18, 732.
18. L. M. Hunnisett, et al., The seventh blind test of crystal structure prediction: structure generation methods, Acta Crystallogr., Sect. B:Struct. Sci., Cryst. Eng. Mater., 2024, 80, 517.
19. X. Liu, J. Zhang and Z. Pei, Machine learning for high-entropy alloys: Progress, challenges and opportunities, Prog. Mater. Sci., 2023, 131, 101018.
20. A. van de Walle and G. Ceder, The effect of lattice vibrations on substitutional alloy thermodynamics, Rev. Mod. Phys., 2002, 74, 11.
21. J. Yang, Mapping temperature-dependent energy-structure-property relationships for solid solutions of inorganic halide perovskites, J. Mater. Chem. C, 2020, 8, 16815.
22. F. Pan, J. Zhai, J. Chen, L. Yang, H. Dong, F. Yuan, Z. Jiang, W. Ren, Z.-G. Ye, G.-X. Zhang and J. Li, Mixed-Halide Perovskite Alloys CsPb(I_1−xBr_x)₃ and CsPb(Br_1−xCl_x)₃: New Insight of Configurational Entropy Effect from First-Principles and Phase Diagrams, Chem. Mater., 2024, 36, 3957.
23. S. Xia, X. Wu, Z. Zhang, Y. Cui and W. Liu, Practical challenges and future perspectives of all-solid-state lithium-metal batteries, Chem, 2019, 5, 753.
24. Q. Liu, Z. Geng, C. Han, Y. Fu, S. Li, Y. B. He, F. Kang and B. Li, Challenges and perspectives of garnet solid electrolytes for all solid-state lithium batteries, J. Power Sources, 2018, 389, 120.
25. C. Wang, K. Fu, S. P. Kammampata, D. W. McOwen, A. J. Samson, L. Zhang, G. T. Hitz, A. M. Nolan, E. D. Wachsman, Y. Mo, V. Thangadurai and L. Hu, Garnet-type solid-state electrolytes: materials, interfaces, and batteries, Chem. Rev., 2020, 120, 4257.
26. F. Han, Y. Zhu, X. He, Y. Mo and C. Wang, Electrochemical stability of Li₁₀GeP₂S₁₂ and Li₇La₃Zr₂O₁₂ solid electrolytes, Adv. Energy Mater., 2016, 6, 1501590.
27. Y. Wang and W. Lai, Phase transition in lithium garnet oxide ionic conductors Li₇La₃Zr₂O₁₂: The role of Ta substitution and H₂O/CO₂ exposure, J. Power Sources, 2015, 275, 612.
28. S. Muy, R. Schlem, Y. Shao-Horn and W. G. Zeier, Phonon-ion interactions: designing ion mobility based on lattice dynamics, Adv. Energy Mater., 2021, 11, 2002787.
29. Y. Chen, E. Rangasamy, C. R. dela Cruz, C. Liang and K. An, A study of suppressed formation of low-conductivity phases in doped Li₇La₃Zr₂O₁₂ garnets by in situ neutron diffraction, J. Mater. Chem. A, 2015, 3, 22868.
30. Z. Yan and Y. Zhu, Impact of Lithium Nonstoichiometry on Ionic Diffusion in Tetragonal Garnet-Type Li₇La₃Zr₂O₁₂, Chem. Mater., 2024, 36, 11551.
31. M. Kubicek, A. Wachter-Welzl, D. Rettenwander, R. Wagner, S. Berendts, R. Uecker, G. Amthauer, H. Hutter and J. Fleig, Oxygen vacancies in fast lithium-ion conducting garnets, Chem. Mater., 2017, 29, 7189.
32. A. Paolella, W. Zhu, G. Bertoni, S. Savoie, Z. Feng, H. Demers, V. Gariepy, G. Girard, E. Rivard, N. Delaporte, A. Guerfi, H. Lorrmann, C. George and K. Zaghib, Discovering the influence of lithium loss on garnet Li₇La₃Zr₂O₁₂ electrolyte phase stability, ACS Appl. Energy Mater., 2020, 3, 3415.
33. D. J. Wales and J. P. Doye, Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms, J. Phys. Chem. A, 1997, 101, 5111.
34. J. Holland, T. Demeyere, A. Bhandari, F. Hanke, V. Milman and C.-K. Skylaris, A Workflow for Identifying Viable Crystal Structures with Partially Occupied Sites Applied to the Solid Electrolyte Cubic Li₇La₃Zr₂O₁₂, J. Phys. Chem. Lett., 2023, 14, 10257.
35. E. P. George, D. Raabe and R. O. Ritchie, High-entropy alloys, Nat. Rev. Mater., 2019, 4, 515.
36. A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater., 2015, 108, 1.
37. S. Adams and R. P. Rao, Ion transport and phase transition in Li_7−xLa₃ (Zr_2−xM_x)O₁₂ (M = Ta⁵⁺, Nb⁵⁺, x = 0,0.25), J. Mater. Chem., 2012, 22, 1426.
38. Y. Wang, A. Huq and W. Lai, Insight into lithium distribution in lithium-stuffed garnet oxides through neutron diffraction and atomistic simulation: Li_7−xLa₃Zr_2−xTa_xO₁₂ (x = 0–2) series, Solid State Ionics, 2014, 255, 39.
39. R. Jalem, M. Rushton, W. Manalastas Jr, M. Nakayama, T. Kasuga, J. A. Kilner and R. W. Grimes, Effects of gallium doping in garnet-type Li₇La₃Zr₂O₁₂ solid electrolytes, Chem. Mater., 2015, 27, 2821.
40. F. A. Garcia Daza, M. R. Bonilla, A. Llordés, J. Carrasco and E. Akhmatskaya, Atomistic insight into ion transport and conductivity in Ga/Al-substituted Li₇La₃Zr₂O₁₂ solid electrolytes, ACS Appl. Mater. Interfaces, 2018, 11, 753.
41. M. Klenk and W. Lai, Local structure and dynamics of lithium garnet ionic conductors: tetragonal and cubic Li₇La₃Zr₂O₇, Phys. Chem. Chem. Phys., 2015, 17, 8758.
42. F. Chen, J. Li, Z. Huang, Y. Yang, Q. Shen and L. Zhang, Origin of the phase transition in lithium garnets, J. Phys. Chem. C, 2018, 122, 1963.
43. M. J. Klenk and W. Lai, Finite-size effects on the molecular dynamics simulation of fast-ion conductors: A case study of lithium garnet oxide Li₇La₃Zr₂O₁₂, Solid State Ionics, 2016, 289, 143.
44. J. Dai, Q. Chen, T. Glossmann and W. Lai, Comparison of interatomic potential models on the molecular dynamics simulation of fast-ion conductors: a case study of a Li garnet oxide Li₇La₃Zr₂O₁₂, Comput. Mater. Sci., 2019, 162, 333.
45. M. R. Bonilla, F. A. G. Daza, J. Carrasco and E. Akhmatskaya, Exploring Li-ion conductivity in cubic, tetragonal and mixed-phase Al-substituted Li₇La₃Zr₂O₁₂ using atomistic simulations and effective medium theory, Acta Mater., 2019, 175, 426.
46. Z. Zhu and Y. Zhu, Developing Classical Interatomic Potentials for Solid Electrolytes, Acc. Mater. Res., 2022, 3, 1101.
47. H. Shiiba, M. Koyama, N. Zettsu and K. Teshima, Li-Ion Conduction Characteristics at Grain Boundaries in Garnet Li_7−xLa₃Zr_2−xNb_xO₁₂ (0 ≤ x ≤ 2), Chem. Mater., 2024, 36, 6370.
48. D. Rettenwander, G. Redhammer, F. Preishuber-Pflügl, L. Cheng, L. Miara, R. Wagner, A. Welzl, E. Suard, M. M. Doeff, M. Wilkening, J. Fleig and G. Amthauer, Structural and electrochemical consequences of Al and Ga cosubstitution in Li₇La₃Zr₂O₁₂ solid electrolytes, Chem. Mater., 2016, 28, 2384.
49. J. Behler, Four generations of high-dimensional neural network potentials, Chem. Rev., 2021, 121, 10037.
50. A. P. Bartók, M. C. Payne, R. Kondor and G. Csányi, Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons, Phys. Rev. Lett., 2010, 104, 136403.
51. A. P. Bartók, R. Kondor and G. Csányi, On representing chemical environments, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 87, 184115.
52. V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti and G. Csányi, Gaussian process regression for materials and molecules, Chem. Rev., 2021, 121, 10073.
53. M. Rupp, A. Tkatchenko, K.-R. Müller and O. A. Von Lilienfeld, Fast and accurate modeling of molecular atomization energies with machine learning, Phys. Rev. Lett., 2012, 108, 058301.
54. Z. Li, J. R. Kermode and A. De Vita, Molecular dynamics with on-the-fly machine learning of quantum-mechanical forces, Phys. Rev. Lett., 2015, 114, 096405.
55. R. Jinnouchi, F. Karsai and G. Kresse, On-the-fly machine learning force field generation: Application to melting points, Phys. Rev. B, 2019, 100, 014105.
56. S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt and K.-R. Müller, Machine learning of accurate energy-conserving molecular force fields, Sci. Adv., 2017, 3, e1603015.
57. J. S. Smith, O. Isayev and A. E. Roitberg, ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost, Chem. Sci., 2017, 8, 3192.
58. S. Chmiela, H. E. Sauceda, K.-R. Müller and A. Tkatchenko, Towards exact molecular dynamics simulations with machine-learned force fields, Nat. Commun., 2018, 9, 3887.
59. K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko and K.-R. Müller, Schnet-a deep learning architecture for molecules and materials, J. Chem. Phys., 2018, 148, 241722.
60. O. T. Unke and M. Meuwly, PhysNet: A neural network for predicting energies, forces, dipole moments, and partial charges, J. Chem. Theory Comput., 2019, 15, 3678.
61. I. Batatia, D. P. Kovacs, G. Simm, C. Ortner and G. Csányi, MACE: Higher order equivariant message passing neural networks for fast and accurate force fields, Advances in Neural Information Processing Systems, 2022, 35, 11423.
62. T. Frank, O. Unke and K.-R. Müller, So3krates: Equivariant attention for interactions on arbitrary length-scales in molecular systems, Advances in Neural Information Processing Systems, 2022, 35, 29400.
63. A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth and B. Kozinsky, Learning local equivariant representations for large-scale atomistic dynamics, Nat. Commun., 2023, 14, 579.
64. B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel and G. Ceder, CHGNet as a pretrained universal neural network potential for charge-informed atomistic modelling, Nat. Mach. Intell., 2023, 5, 1031.
65. I. Batatia, et al., A foundation model for atomistic materials chemistry, arXiv, 2023, preprint, arXiv:2401.00096.
66. A. C. Grieder, K. Kim, L. F. Wan, J. Chapman, B. C. Wood and N. Adelstein, Effects of Nonequilibrium Atomic Structure on Ionic Diffusivity in LLZO: A Classical and Machine Learning Molecular Dynamics Study, J. Phys. Chem. C, 2024, 128, 8560.
67. K. Kim, A. Dive, A. Grieder, N. Adelstein, S. Kang, L. F. Wan and B. C. Wood, Flexible machine-learning interatomic potential for simulating structural disordering behavior of Li₇La₃Zr₂O₁₂ solid electrolytes, J. Chem. Phys., 2022, 156, 221101.
68. M. Ataya, E. McCalla and R. Z. Khaliullin, Machine Learning for High-Throughput Configuration Sampling of Li-La-Ti-O Disordered Solid-State Electrolyte, J. Phys. Chem. C, 2024, 128, 14149.
69. A. Duval, S. V. Mathis, C. K. Joshi, V. Schmidt, S. Miret, F. D. Malliaros, T. Cohen, P. Lio, Y. Bengio and M. Bronstein, A Hitchhiker's Guide to Geometric GNNs for 3D Atomic Systems, arXiv, 2023, preprint, arXiv:2312.07511.
70. J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals and G. E. Dahl, International Conference in Machine Learning, 2017, 1263.
71. S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nat. Commun., 2022, 13, 2453.
72. J. Wang, Y. Wang, H. Zhang, Z. Yang, Z. Liang, J. Shi, H.-T. Wang, D. Xing and J. Sun, E(n)-Equivariant cartesian tensor message passing interatomic potential, Nat. Commun., 2024, 15, 7607.
73. X. Fu, Z. Wu, W. Wang, T. Xie, S. Keten, R. Gomez-Bombarelli and T. Jaakkola, Forces are not enough: Benchmark and critical evaluation for machine learning force fields with molecular simulations, arXiv, 2022, preprint, arXiv:2210.07237.
74. J. T. Frank, O. T. Unke, K.-R. Müller and S. Chmiela, A Euclidean transformer for fast and stable machine learned force fields, Nat. Commun., 2024, 15, 6539.
75. A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser and I. Polosukhin, Attention Is All You Need, arXiv, 2023, preprint, arXiv:1706.03762.
76. P. Korotaev, M. Belov and A. Yanilkin, Reproducibility of vibrational free energy by different methods, Comput. Mater. Sci., 2018, 150, 47.
77. D.-B. Zhang, T. Sun and R. M. Wentzcovitch, Phonon quasiparticles and anharmonic free energy in complex systems, Phys. Rev. Lett., 2014, 112, 058501.
78. V. S. Kandagal, M. D. Bharadwaj and U. V. Waghmare, Theoretical prediction of a highly conducting solid electrolyte for sodium batteries: Na₁₀GeP₂S₁₂, J. Mater. Chem. A, 2015, 3, 12992.
79. A. H. Larsen, et al., The atomic simulation environment a Python library for working with atoms, J. Phys.:Condens. Matter, 2017, 29, 273002.
80. B. Ahammed and E. Ertekin, Configurational Disorder, Strong Anharmonicity, and Coupled Host Dynamics Lead to Superionic Transport in Li₃YCl₆ (LYC), Adv. Mater., 2024, 36, 2310537.
81. M. Baggioli and A. Zaccone, Universal origin of boson peak vibrational anomalies in ordered crystals and in amorphous materials, Phys. Rev. Lett., 2019, 122, 145501.
82. F. Knoop, T. A. Purcell, M. Scheffler and C. Carbogno, Anharmonicity measure for materials, Phys. Rev. Mater., 2020, 4, 083809.
83. J. Yang, J. Fan and S. Li, Mapping the Room-Temperature Dynamic Stabilities of Inorganic Halide Double Perovskites, Chem. Mater., 2022, 34, 9072.
84. T. D. Humphries, J. Yang, R. A. Mole, M. Paskevicius, J. E. Bird, M. R. Rowles, M. S. Tortoza, M. V. Sofianos, D. Yu and C. E. Buckley, Fluorine substitution in magnesium hydride as a tool for thermodynamic control, J. Phys. Chem. C, 2020, 124, 9109.
85. F. Noé, S. Olsson, J. Köhler and H. Wu, Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning, Science, 2019, 365, eaaw1147.
86. S. Mehdi, Z. Smith, L. Herron, Z. Zou and P. Tiwary, Enhanced sampling with machine learning, Ann. Rev. Phys. Chem., 2024, 75, 347.
87. L. Midgley, V. Stimper, J. Antorán, E. Mathieu, B. Schölkopf and J. M. Hernández-Lobato, SE (3) equivariant augmented coupling flows, Advances in Neural Information Processing Systems, 2024, 36, 79200.
88. F. H. Stillinger and T. A. Weber, Packing structures and transitions in liquids and solids, Science, 1984, 225, 983.
89. D. J. Wales, Coexistence in small inert gas clusters, Mol. Phys., 1993, 78, 151–171.
90. B. Strodel and D. J. Wales, Free energy surfaces from an extended harmonic superposition approach and kinetics for alanine dipeptide, Chem. Phys. Lett., 2008, 466, 105.
91. V. A. Sharapov, D. Meluzzi and V. A. Mandelshtam, Low-temperature structural transitions: Circumventing the broken-ergodicity problem, Phys. Rev. Lett., 2007, 98, 105701.
92. Y. Zhu, J. Zhang, W. Li, Y. Zeng, W. Wang, Z. Yin, B. Hao, Q. Meng, Y. Xue, J. Yang and S. Li, Enhanced Li⁺ conductivity of Li₇La₃Zr₂O₁₂ by increasing lattice entropy and atomic redistribution via Spark Plasma Sintering, J. Alloys Compd., 2023, 967, 171666.
93. Y. Zhu, J. Zhang, W. Li, Y. Xue, J. Yang and S. Li, Realization of superior ionic conductivity by manipulating the atomic rearrangement in Al-doped Li₇La₃Zr₂O₁₂, Ceram. Int., 2023, 49, 10462.
94. E. Anderson, E. Zolfaghar, A. Jonderian, R. Z. Khaliullin and E. McCalla, Comprehensive Dopant Screening in Li₇La₃Zr₂O₁₂ Garnet Solid Electrolyte, Adv. Energy Mater., 2024, 14, 2304025.
95. J. A. Harrison, J. D. Schall, S. Maskey, P. T. Mikulski, M. T. Knippenberg and B. H. Morrow, Review of force fields and intermolecular potentials used in atomistic computational materials research, Appl. Phys. Rev., 2018, 5, 031104.
96. D. Minh, H. X. Wang, Y. F. Li and T. N. Nguyen, Explainable artificial intelligence: a comprehensive review, Artif. Intell. Rev., 2022, 55, 1.
97. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169.
98. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758.
99. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett., 2008, 100, 136406.
100. S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson and G. Ceder, Python Materials Genomics (pymatgen): A robust, open-source Python library for materials analysis, Comput. Mater. Sci., 2013, 68, 314.
101. F. Musil, S. De, J. Yang, J. E. Campbell, G. M. Day and M. Ceriotti, Machine learning for the structure-energy-property landscapes of molecular crystals, Chem. Sci., 2018, 9, 1289–1300.
102. D. P. Kingma, Adam: A method for stochastic optimization, arXiv, 2014, preprint, arXiv:1412.6980.
103. D. C. Liu and J. Nocedal, On the limited memory BFGS method for large scale optimization, Math. Prog., 1989, 45, 503.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp00138b

‡ By convention, the letters represent the Wyckoff positions, and the numbers signify the total number of equivalent Wyckoff positions, in the Ia3̄d space group for c-LLZO.

§ In the numerical implementation of eqn (2), H⁽ⁱ⁾ is measured with respect to the energy zero point, which, in this case, is set at min{H⁽ⁱ⁾}, namely, the global energy minimum for the set of LLZO under consideration. The same principle also applies to eqn (4).

---