Accurate Surface and Finite-Temperature Bulk Properties of Lithium Metal at Large Scales Using Machine Learning Interaction Potentials

The properties of lithium metal are key parameters in the design of lithium-ion and lithium-metal batteries. They are difficult to probe experimentally due to the high reactivity and low melting point of lithium as well as the microscopic scales at which lithium exists in batteries where it is found to have enhanced strength, with implications for dendrite suppression strategies. Computationally, there is a lack of empirical potentials that are consistently quantitatively accurate across all properties, and ab initio calculations are too costly. In this work, we train a machine learning interaction potential on density functional theory (DFT) data to state-of-the-art accuracy in reproducing experimental and ab initio results across a wide range of simulations at large length and time scales. We accurately predict thermodynamic properties, phonon spectra, temperature dependence of elastic constants, and various surface properties inaccessible using DFT. We establish that there exists a weak Bell–Evans–Polanyi relation correlating the self-adsorption energy and the minimum surface diffusion barrier for high Miller index facets.


Introduction
Lithium metal batteries provide a promising pathway to achieving high capacity energy storage devices.[3] A number of approaches have been proposed to address the issue of dendrite formation.
One approach is suppressing instability through the introduction of a solid electrolyte in contact with lithium metal.Monroe and Newman proposed the use of a solid polymer electrolyte with a shear modulus larger than that of lithium, 4 which was extended by Ahmad and Viswanathan showing that an alternate approach could be to use a lithium-dense solid electrolyte whose modulus is smaller than that of lithium. 5A second approach increases morphological stability through rapid surface diffusion, quantified by surface diffusion barriers. 6,7All of these approaches critically hinge upon accurate determination of the properties of lithium metal.The first approach requires knowing the room temperature mechanical properties of lithium, while the second approach requires a detailed understanding of surface diffusion barriers across high Miller index facets formed during morpohological instability.
Determining these properties experimentally is challenging due to the high reactivity of lithium 8 thus computational methods provide a more practical approach.
In principle, many material properties can be calculated from first principles using high fidelity, atomistic methods such as Ab-Initio Molecular Dynamics (AIMD).In practice however, the computational cost of these methods is often prohibitive due to poor scaling with the number of electrons (∼O(N 3 )) and the long simulation time per timestep. 9As an alternative, approximate empirical potentials have been commonly used for Molecular Dynamics (MD) simulations to achieve necessary time and length scales to collect good statistics but often at the cost of considerable loss of accuracy leading to only qualitative or often wrong results for some properties.[12] MLIPs have been used to study supercritical phenomena in hydrogen, 13 defects in various metals, 14 have been benchmarked for transition metals 15 and many other examples.For lithium, the SNAP potential 16 has been developed for purposes of benchmarking MLIPs and is therefore limited in its applicability as we show in this work.Jiao et al. generated a Deep Potential 17 and simulated the self-deposition of Li and the different morphologies that could arise in deposition processes. 18Their potential however was not accurate in predicting stresses and elastic constants, which we improve upon in our own Deep Potential.
In this work, we generate data to train two general purpose MLIPs for pure lithium metal based on NequIP 11 and Deep Potential.The MLIPs, particularly NequIP, reproduce DFT and experimental results remarkably well over a wide range of structures including, bulk, surfaces, defects and liquids, all in one potential, consistently outperforming empirical potentials and existing MLIPs.We therefore more accurately calculate elastic and surface properties important to the design of lithium metal batteries and discuss the implications of our results.

Results
In this section, we first demonstrate the accuracy of the trained MLIPs by comparing the predictions directly to results from DFT and other potentials in the literature.Second, we present results from a number of simulations that are typically too expensive or otherwise impossible to do with DFT and compare them to experiment.Finally we predict key properties in the design of lithium metal batteries using simulations that have, up to this point, been inaccessible using DFT using the NequIP model architecture.
In addition to DFT, where possible, we compare our results with the popular MEAM empirical potential developed by Kim et al. 19 and SNAP MLIP by Zuo et al. 16 The MEAM potential has been used to predict, for example, thermal behavior in lithium metal electrodes 20 and the lifetime of glassy lithium nuclei under fast charging conditions. 21SNAP is a MLIP that was designed for benchmarking a number of MLIP architectures and not necessarily built for production simulations.It is worth noting that the differences in predictions with our MLIPs can be attributed not only to the quality of the fit but also to the difference in the datasets used to fit the potentials and therefore the potentials need to be compared with experiment to be assessed rigorously.We perform many of these assessments in the following sections.
We show results for two NequIP potentials trained with different levels of precision.
NequIP32 and NequIP64 are trained to single (float32) and double (float64) precision respectively as Batatia et al. found that float32 can be insufficiently precise to represent total energy differences. 22In our experiments, we found that the lower precision of NequIP32 led to a more coarse discretization of the potential energy surface resulting in errors of ∼10meV for a modest number of ∼100 atoms when calculating energy differences.The NequIP32 potential however uses less memory and is faster hence it would be advantageous to use a lower precision, particularly for MD simulations which depend on forces and are not affected by small errors in energy and do not have changing numbers of atoms and we have found to be stable.Properties like the vacancy formation energies however are affected by this precision hence we show both results.We found that the numerical instabilities were improved by performing operations in the last layers of the model architecture in double precision, however we did not implement this solution for the potentials in this work.We present results for both NequIP32 and NequIP64 to highlight any differences.All MD simulations however were performed at single precision to save computational effort.We also show results for our developed Deep Potential (DP).All calculations for DP were performed at double precision.Predictions for a number of different properties accessible using standard DFT methods are shown in Table 1 to benchmark the MLIPs.Details on how the calculations were performed are given in the Supporting Information (SI Appendix).

DFT Benchmarks
The conventional metric to assess the accuracy of an MLIP is the Root Mean Square Error (RMSE) on the prediction of the label y which can be energy, forces and stresses over a test set.The RMSE is defined as where y i is the ground truth label and ŷi is the predicted label of sample i in a test set with N samples.The RMSEs shown in Table 1 demonstrate very good accuracy on the order of ∼1meV/atom for energies, comparable to the fluctuations in the energy with converged k points and typically chosen DFT convergence criteria. 24The forces and stresses are also well converged to ∼10meV/ Å and ∼1meV/ Å3 respectively such that an atomic/cell parameter displacement of ∼ 0.01 Å gives an energy difference of ∼1meV/atom.It is worth noting that both NequIP potentials perform an order of magnitude better in stress error (0.06eV/ Å3 ) compared to DP (0.22eV/ Å3 ) and are slightly better than DP in most other metrics.We also note that we trained the NequIP architecture on the data that was used to train the SNAP potential and got better test errors on the SNAP test set with little hyperparameter tuning.
NequIP RMSEs were 1.3meV/atom and 15mev/ Å vs quoted SNAP RMSEs of 1.4meV/atom and and 40meV/ Å for SNAP for energy and force RMSE respectively. 16We thereby expect that the NequIP architecture is generally more accurate for the same dataset, consistent with previous results. 11e RMSE is a good metric for benchmarking MLIP architectures against each other but does not imply an accurate potential in predicting experimental properties since it strongly depends on an inevitably biased test set.We therefore calculate a number of other properties to benchmark the MLIPs which reflect the quality of the distribution of data in the training set as well as the appropriateness of the choice of ground-truth data, particularly the choice of DFT exchange-correlation functional.Table 1 shows how all the properties calculated here are in excellent agreement with most being within 5% of the DFT predictions for the DP and NequIP potentials.An exception is the anisotropy which is derived from the elastic constants as 2C 44 /(C 11 − C 12 ) and hence propagates error from the elastic constants.

Phonon Spectrum
We calculate the phonon spectrum of BCC lithium using the phonons package as implemented in the Atomic Simulation Environment (ASE) in Fig. 1a.The finite displacement method was used with a displacement of 0.01 Å for a supercell of size 5x5x5.Note that while this supercell is larger than the cutoff in the MLIPs of 6 Å, the NequIP potentials can still predict slight changes in force due to their message passing architecture with multiple layers, resulting in a much larger effective cutoff radius capable of capturing some long range effects.for some parts of the Brillouin zone. 26The temperature effects can be taken into account by sampling the dynamical matrix at finite temperature but this is beyond the scope of this work.The structural features of the bands however are clearly well captured by NequIP, DP and SNAP which produces slightly different results, including the broadening/splitting of the phonons between Γ-H, typical of anisotropic cubic materials like lithium. 26MEAM however completely fails to capture the phonon spectrum as it was not specifically trained on phonon data due to having a short cutoff that does not capture long range interactions important at small wave vectors.The Modified Analytic Embedde Atom Method (MAEAM) has been found to predict a reasonable spectrum. 27

Surface Energies and Wulff Construction
As shown in Fig. 1b, in addition to bulk properties, the NequIP and DP potentials accurately describe the surface energies of lithium.Surface and adsorption energies are often quoted up to ∼10meV/ Å2 precision due to finite size effects and propagation of errors.The NequIP and DP potentials are well within 1meV/ Å2 error for all the miller indices despite being only explicitly trained on the (100), ( 110) and (111) planes as starting seeds.This demonstrates the excellent generalization to higher miller index surfaces which we take advantage of in the calculation of surface properties in Section .The predicted surface energies for low miller indices listed in Table 1 agree very well with both our DFT and DFT results in other works whereas the MEAM and SNAP predictions give different values and slightly different rankings for the surface energies as a function of miller index suggesting they are less reliable for surface properties.
We also calculate the Wulff construction which determines the equilibrium shape of a droplet or crystal suspended in a medium which in this case is vacuum. 28Depending on the surface energies and geometry, the shape of the droplet determines what fraction of the total area of the droplet is contributed by each facet.This allows us to estimate the importance of considering a particular facet in simulations at equilibrium.In the presence of an electrolyte, the Wulff construction might change due to the sensitivity to the small energy differences between higher miller indices and in practice the surrounding medium. 23We follow the procedure used by Tran et al. 23 implemented in the Python Materials Genomics package (pymatgen) to estimate the area fractions for miller indices up to a value of 3 in Fig. 1c.
The NequIP results agree well with the DFT prediction that (110) and (100) planes followed by ( 111) and (211) to a lesser extent dominate the Wulff construction.DP has slightly different results with (211) more dominant over (111) while MEAM and SNAP have much more significant differences as expected due to the poor prediciton of surface energies.The domination of (110) in the Wulff construction despite (100) having the lowest surface energy is consistent with experimental results. 8e excellent reproduction of DFT calculated properties gives confidence that the MLIPs, particularly NequIP are truly reproducing the DFT result for surface and bulk properties with small errors at the level of typical DFT precision.We therefore assume subsequent errors are from the quality of the dataset, most likely the choice of exchange-correlation functional when comparing with experiment.In the rest of the work, we demonstrate the superior accuracy in reproducing experimental results of the NequIP potentials over a range of properties that are difficult to predict using DFT.

Temperature Dependence of Elastic Properties of Lithium
We estimate a number of mechanical properties using MD implemented in the Large-scale Figure 2: Various bulk mechanical properties of lithium as a function of temperature.The NequIP and DP potentials in this work perform well across all bulk properties for which there is significant experimental data while MEAM performs reasonably well but slightly worse and SNAP performs poorly.We also compare the predictions to the Quasiharmonic Approximation (QHA) by 29  All the MLIPs and MEAM perform well with less than 1.5% error for the lattice constant with a slight underestimation for the MLIPs.MEAM overestimates the potential as the 0K parameters are fit to match room temperature result of 3.50 Å. 19 The temperature range considered was chosen because there exists a Martensitic transition into a FCC structure at 78K and the melting point of lithium is 450K. 26NequIP, MEAM and DP predict a strongly linear dependence of the lattice constant with temperature whereas SNAP predicts a slightly increasing thermal expansion coefficient corresponding to an increasing slope of the lattice constant with temperature.There is not enough experimental data above room temperature to determine which behavior is most likely.Such experiments may be challenging due to plasticity as they would be performed at a high homologous temperature, greater than 0.7 times the melting point for lithium.
The elastic response of single crystal BCC lithium, particularly at microscopic scales is key to the design of LMSSBs.The bulk and shear modulus are key parameters in the model of Monroe and Newman as well as Ahmad and Viswanathan used to predict stability against the formation of dendrites. 4,5Due to lithium's low melting point, it is a soft material whose mechanical response can have unique properties near room temperature.The interplay between elastic and plastic regimes has been a topic of study. 29,34,35Here, we calculate the elastic constants of single crystal BCC lithium as a function of temperature.Xu et al.
measured the elastic constants of lithium nanopillars and proceeded to calculate bulk elastic constants using a Quasiharmonic approximation (QHA) within DFT. 29 They found that the QHA performed poorly hence the need for simulations at the fidelity of AIMD.
We perform the calculations for the elastic constants C 11 , C 12 and C 44 by fitting stressstrain curves after applying two different types of strain, an orthorhombic and monoclinic strain and perform 1ns long NVT simulations following the prescription by Zhang et al. 36 to save computational cost.All other elastic constants for cubic crystals can be derived from these three using well known formulae 37 implemented in pymatgen. 38e predictions of the components of the C 11 , C 12 and C 44 , the Universal Anisotropy, Voigt-Reuss-Hill averaged bulk and shear moduli (K VRH and G VRH respectively) as well as the Young's modulus and Poisson Ratio are plotted as a function of temperature in Fig. 2. The Voigt-Reuss-Hill (VRH) average is a reliable method for predicting polycrystalline elastic moduli given the relevant single crystal elastic constants while being very easy to compute.37 As shown in Fig. 2, the MLIPs in this work reproduce the experimental results well, significantly outperforming the QHA, SNAP and MEAM in reproducing experimental results for single crystal lithium and in the prediction of VRH averaged quantities.
The QHA is the only model that fails to predict the C 44 qualitatively accurately underestimating the dependence of C 44 as a function of temperature, potentially due to the assumed independence of vibrational modes in each spatial dimension.SNAP performs very poorly on C 11 , C 12 , predicting that they increase before decreasing with temperature, inconsistent with all experimental and previous computational results.The properties derived from these constants accumulate errors therefore SNAP is not suitable for calculations that rely on stress prediction.MEAM performs generally well on all predictions as it was specifically fit to elastic constants but predicts an increasing anisotropy with temperature different from all the other potentials and experiment.
DP and NequIP32 are in excellent agreement with experimental results, consistently within 10% or less of the experimental results for C 11 , C 12 and C 44 and with matching qualitative behavior.NequIP32 performs slightly better than DP, likely due to the more accurate stress predictions, but both potentials are clearly suitable for prediction of bulk elastic properties.The only exception is the Anisotropy which is overestimated by DP and NequIP32 at low temperatures and shows much larger decrease with increasing temperature than experimental predictions.
Overall, the DP and NequIP32 potentials in this work are the most accurate potentials with which to calculate bulk and elastic phenomena for BCC lithium.

Adsorption and Surface Diffusion Barriers
Another set of key properties for modeling lithium in the context of batteries involves surface properties such as adsorption and diffusion across the surface.There are a number of studies that perform these simulations using DFT but they are often forced to limit the surface area of the slabs, the miller indices considered and the number of layers in their slab models. 6,39,40 our experiments, we found that it was impossible to converge the number of layers in the calculation of adsorption energy to within ∼1meV/atom, even for the lowest miller indices, with less than 10 layers.Sun et al. noted that the poor convergence of surface energies was likely due to inconsistencies in the Brillouin zone sampling of the clean slabs, bulk references and adsorbed slabs. 41They introduced a workaround where a slab-consistent bulk reference is used for faster convergence with fewer layers for surface energy.However, it is not clear how to create a bulk reference consistent with an adsorbed surface hence a need for the surface energy to converge even with a single atom bulk reference for adsorption calculations.Since the MLIP is not limited as much as DFT in the number of layers, we can adequately converge the number of layers, typically on the order of tens of layers even with large surface areas.
For DFT calculations of adsorption energies and surface diffusion barriers, the workaround is to fix a few layers, usually 2 at the bottom of a slab with 4-6 layers in total, artificially imposing a bulk environment close to the surface. 24This limits the possible facets which can be considered to the very lowest few miller indices, typically up to 211.This is hampered as well by the failure of convergence with so few layers.With the MLIPs, the size of the slab model is less restricted thus we construct slabs of at least 14 Å in thickness with the first 6 Å of layers fixed consistent with the exposed facet onto which an adsorbate is placed in all our calculations using pymatgen.
For all results in this section, we use the most accurate NequIP64 potential.The DP potential also performs well on all the properties but to a lesser extent.
We plot the Surface Potential Energy Surface (SPES) for various miller indices in Fig. 3a.The SPES is calculated by adsorbing a lithium atom on the relevant facet at various positions in the surface unit cell of a 4x4 surface supercell and allowing the adatom to relax only in the perpendicular direction (z direction) to the surface while allowing all other atoms (except the 6 fixed layers at the bottom of the slab) to relax in any direction.The adsorption energy (E ads ) is calculated using an unconventional method that uses a bulk reference for the adatom as The bulk reference is chosen as it is more reliably predicted by all the potentials since no single atom data was added to the training datasets and we are only interested in relative adsorption energies and activation barriers.Relative Adsorption Energy is used in the SPES i.e. the energy above the lowest adsorption energy for that particular facet.The SPES shows the diverse geometries and distributions of adsorption sites that can arise with exposure of different facets and can potentially be used for thermodynamic lattice models of adsorption and surface diffusion properties.SPESs for higher miller indices can be found in the Supporting Information (SI Appendix).
Using the same slab models as those used for adsorption energies, we also calculate the surface diffusion barriers as an adatom diffuses from one surface unit cell to the next using the Nudged Elastic Band (NEB) method with a spring constant of 0.1eV/ Å and 7 images as implemented in ASE. 42The paths with the lowest diffusion barriers identified by the NEB method are shown on the SPES in Fig. 3a.We find that for higher miller indices, the lithium atoms are more likely to diffuse via one-dimensional channels formed by step edges as there are higher coordination numbers in this environment and hence lower adsorption energies in agreement with Gaissmaier et al.Only (100) and ( 110) and (111) have equally high diffusion barriers in two dimensions. 7 Fig. 3b, we plot the adsorption energy versus the lowest surface diffusion barrier for that facet and show that there exists a linear correlation between the adsorption energy and the surface diffusion barrier.This is an example of a Bell-Evans-Polanyi (BEP) principle, a general framework that notes that the difference in activation energy between two reactions of the same family is roughly proportional to the difference of their enthalpy of reaction. 43,44 the context of adsorption energies and surface diffusion barriers, Pande et al. found a BEP relation relating the adsorption energies on low miller index facets for different metals and alloys and the surface diffusion barriers on those surfaces. 40They argued that the BEP relation can be used as a powerful screening tool that avoids expensive NEB calculations.
We show that there also exists a BEP relation for different miller indices of the same material which can similarly be used for screening purposes in cases where the NEB method is too costly.We note that the (111) surface is an outlier to the BEP relation, likely due to it's very high surface energy as it is not a close-packed surface for the BCC crystal structure and all the surface atoms are under-coordinated.We also predict an example of an Erlich-Schwöebel barrier, a descriptor used in determining the likelihood of an adatom to diffuse up a step-edge. 45The activation energy for this process is increased as the adatom goes from a region of high coordination along the step-edge to lower coordination on top of the step.The Erlich-Schwöebel activation barrier (E S ) is defined as where E ESB is the increased barrier due to having to diffuse up a step-edge and E T is the activation energy for diffusion without any steps. 7We consider two mechanisms for step-up diffusion on the most stable (100) surface and show in Fig. 3c that NequIP64 reproduces the results by Gaissmaier et al. to within 0.01eV for the activation energy for both mechanisms. 7schematic to explain the two mechanisms is shown in Fig. 3c.The first mechanism is when the adatom directly diffuses (Diff.)over the step-edge (1→3) and the second is by exchanging (Ex.) with an atom in the step (1→2, 2→3).
Overall the NequIP64 potential has shown very good accuracy in reproducing and predicting surface properties and has allowed us to predict properties that are computationally infeasible using DFT.

Discussion: Implications for Li-metal Battery Design
The calculations enabled by the lithium MLIPs derived in this work now allow refining criterion for morphological stability associated with dendrite formation described above.Monroe and Newman, subsequently, Ahmad and Viswanathan showed that room temperature shear modulus, anisotropy of lithium are critical factors in determining the stability criterion for dendrite suppression and thereby the required shear modulus of the solid electrolyte. 4,5In this work, we show that the temperature-dependent response of the mechanical properties is much stronger than that previously predicted using the QHA.The room-temperature modulus, as predicted by NequIP32 is around 3.5 GPa, a reduction of about ∼30% from that at 0K using the QHA, drastically modifying the stability criterion, thereby modifying the required shear modulus of the solid electrolyte needed to suppress dendrite formation by a similar ∼30%.surface or many orders of magnitude more, even along preferential one-dimensional channels observed in the SPES.
Another implication is that glassy phases which have been shown to improve bulk conductivity, 21 are likely to possess lower surface diffusion limiting the charging rate when glassy phases of lithium are formed.Surface diffusion in glassy phases is further hampered by the lack of long range order as lithium is deposited which would randomly orient the preferential one-dimensional channels with high adsorption energy thereby trapping atoms in deep potential wells.We therefore expect that the combined effect is the proliferation of local instabilities under fast charging conditions before the surface equilibrates to (100) and (110) facets with faster diffusion in two dimensions which would reduce dendrite growth.
We believe that the lithium MLIP and training dataset developed here will enable the calculation of void formation, creep behavior and various other meso-scale properties previously not tractable using atomistic simulations with significantly improved accuracy thereby allowing refined material properties under different conditions.

Methods
In this section, we describe the model architectures, data generation and simulation details.

Machine Learning Interaction Potentials
MLIPs are parametric models that can be trained on a set of atomic configurations given in terms of coordinates {R} and labelled with the corresponding energies, forces and stresses for that configuration from a high fidelity source such as DFT.
Typically, the total energy of an atomic configuration is modeled as a sum of atomic energies as in () i.e.

E({R})
The i-th atom contributes an energy (ε i ) that is learned from the geometry centered on atom i.The forces on each atom and the virial stresses can then be calculated as the derivatives of the total energy with respect to the atomic positions and the cell dimensions respectively using automatic differentiation.It is essential, that the predicted energies, forces and stress are appropriately equivariant with respect to translations, rotations and permutations of atoms of the same species.The key difference between different potentials is how they implement this equivariance as described for Deep Potential and NequIP.

Deep Potential -Deep Potential provides a trainable, fully local framework for training
MLIPs that can predict energy, forces and stresses of a particular configuration of atoms.In this work, we choose the se e3 descriptor constructed from radial and angular information of atomic environments.This descriptor is invariant to the symmetries and not trainable.
The neural network architecture is that of Residual Neural Network. 17 NequIP -The Neural Equivariant Interatomic Potential (NequIP) is an E(3)-equivariant message passing interatomic potential that was shown to demonstrate state-of-the-art accuracy, sample efficiency, and transferability on a variety of materials systems at the time of writing. 11,22,25While conventional interatomic potentials operate on invariant descriptors of the materials systems, such as distances and angles, NequIP directly operates on relative interatomic positions ⃗ r ij represented as a graph and leverages latent features comprised of not only scalar, but also vector and higher-order tensor features.

DFT Data
Bulk: Crystalline/liquid structures sampled at different temperature and pressure

Identify relevant atomic environments as seeds
Surface: Slabs of varying width and coverage, crystalline and molten.
Defected: Interstitial, Vacancies, Grain Boundaries A schematic of the procedure used to generate data is shown in Fig. 4. All data used to train the MLIPs was generated using Density Functional Theory with the same parameters across the entire dataset.The parameters chosen were such that the brillouin zone sampling density, plane wave and density cutoffs were converged to ¡1meV/atom to ensure consistency and give a concrete estimate of statistical noise in the data.DFT calculations were done using Quantum Espresso 47 within the Generalized Gradient Approximation using the Perdew-Burke-Eizenhoff exchange correlation functional 48 and the Projector Augemented Wave approach 49 with a plane wave cutoff energy of 1360eV.We used the pseudopotential Li.pbe-s-kjpaw psl.1.0.0.UPF from http://www.quantum-espresso.org.A uniform Brillouin Zone spacing of 0.02 Å−1 with a Monkhorst-Pack 50 sampling procedure was used.To help with convergence of the the Fermi surface, Methfessel-Paxton 51 smearing using a smearing width of 0.27eV was chosen.
The data used to train the potentials was generated using an active learning approach as implemented in the DPGen package.

Figure 1 :
Figure 1: Various properties of BCC lithium calculated and compared with DFT results.a) Phonon spectrum calculated with the various potentials, NequIP and DP from this work as well as SNAP perform well but MEAM completely fails to predict a reasonable spectrum.b) Surface energies for various facets of BCC lithium.NequIP and DP potentials reproduce DFT results very well while MEAM and SNAP have large errors relative to the DFT prediction.c) The Wulff construction for BCC lithium showing that BCC lithium in vacuum is dominated by the (100), (110), (111) and (211) facets according to DFT.NequIP and DP potentials in this work reproduce that result but SNAP and MEAM predict a significant contribution from other facets and much less from (111).
Atomic/Molecular Massively Parallel Simulator (LAMMPS) using the potentials considered in this work.The results are used to predict key parameters in the design of Lithium Metal Solid State Batteries (LMSSBs) such as the bulk and shear modulus as a function of temperature.The temperature dependence of the lattice constant of BCC lithium is shown in Fig.2a.Starting from the conventional equilibrium BCC structure in a 6x6x6 supercell at 0K for each potential, an NPT simulation is used to raise the temperature at a heating rate of 0.01K/timestep and then run at equilibrium for 100,000 timesteps with 1fs/timestep at Figure2: Various bulk mechanical properties of lithium as a function of temperature.The NequIP and DP potentials in this work perform well across all bulk properties for which there is significant experimental data while MEAM performs reasonably well but slightly worse and SNAP performs poorly.We also compare the predictions to the Quasiharmonic Approximation (QHA) by29 and demonstrate the limitation of the QHA approach in finite temperature elastic constant prediction, especially for C 44 .Error bars in the lattice constant (a) are the standard deviation due to fluctuations in the volume of the NPT simulation.Error bars in the elastic constants (b-d) are standard errors from the fitting of stress-strain curves and errorbars in the anisotropy (e) are propagated from errors in the elastic constants.

Figure 3 :
Figure 3: Various surface properties calculated using the NequIP64 potential.a) Demonstration of a BEP correlation between the adsorption energy and the minimum diffusion barrier of each facet.The (111) surface was not included in the fitting of the dotted line.b) Calculation of an Erlich-Schwöebel barrier matching results in the literature.c) A variety of Surface Potential Energy Surfaces colored by the Relative Adsorption Energy show the the different geometries of adsorption sites and the minimum diffusion barrier paths adatoms can take when diffusing from one surface unit cell to the next.
Jäckle et al. and Gaissmaier et al.only probed surface diffusion barriers for low Millerindex facets, tractable using DFT calculations.6,7Morphological protrusions can often locally have high Miller index domains with a high degree of under-coordination such as the curved surface of a roughly cylindrical dendrite.Using the lithium MLIP developed here, we show the existence of a BEP relation, indicating that high Miller indices which typically have higher lithium adsorption energy due to under-coordination have much larger surface diffusion barriers.For the (110) and (211) facets which appear in the Wulff construction, the diffusion rates ν ∝ exp (−∆E/k B T ) at 300K are 1-2 times slower than on the (100) facet with the lowest barrier.The exception is the (111) surface which has an anomalously large surface diffusion barrier with respect to the BEP relation.Other facets which do not appear in the Wulff construction all have higher diffusion barriers with ν >4 times that of the (100)

Figure 4 :
Figure 4: Schematic showing the process by which data was generated and potentials were validated.
52 A variety of seed structures including bulk crystals, monovacant, mono-interstitial, clean surfaces and surfaces with varying coverage up to a maximum of 128 atoms.The seed structures were used as starting structures for isobaric, isothermal (NPT) MD simulations.The MD simulations were performed sequentially on a grid with temperatures of 50K, 300K, 450K and 900K i.e. double the melting temperature and pressures of 0GPa, 0.1GPa, 1GPa and 10GPa.With each active learning step in the sequence, data from previous steps was incorporated into subsequent steps to improve accuracy and flag new structures with high uncertainty for labeling.The final dataset used to train production potentials was split into 4548 structures in the training set and 505 configurations in the test set.Convergence of the dataset and potentials was determined using RMSEs and accurate property prediction across phenomena in all the different classes of configurations considered relevant for a lithium potential of the entire phase space.

Table 1 :
Various properties of BCC lithium predicted using the MLIPs and compared to DFT results in the literature, DFT in this work and the existing MEAM and SNAP potentials.Percentage errors relative to the DFT prediction are shown in square brackets.Except for the Anisotropy, all the errors are within 5% for the NequIP64 potential.NequIP32 predicts has a large vacancy formation energy error due to the lower precision used to calculate small energy differences.SNAP and MEAM produce different results as they were trained on different data and parameterized differently than in this work.