An Electrostatic Spectral Neighbor Analysis Potential (eSNAP) for Lithium Nitride

Machine-learned interatomic potentials based on local environment descriptors represent a transformative leap over traditional potentials based on rigid functional forms in terms of prediction accuracy. However, a challenge in their application to ionic systems is the treatment of long-ranged electrostatics. Here, we present a highly accurate electrostatic Spectral Neighbor Analysis Potential (eSNAP) for ionic $\alpha$-\ce{Li3N}, a prototypical lithium superionic conductor of interest as a solid electrolyte or coating for rechargeable lithium-ion batteries. We show that the optimized eSNAP model substantially outperforms traditional Coulomb-Buckingham potential in the prediction of energies and forces, as well as various properties, such as lattice constants, elastic constants and phonon dispersion curves. We also demonstrate the application of eSNAP in long-time, large-scale Li diffusion studies in \ce{Li3N}, providing atomistic insights into measures of concerted ionic motion (e.g., the Haven ratio) and grain boundary diffusion. This work aims at providing an approach to developing quantum-accurate force fields for multi-component ionic systems under the SNAP formalism, enabling large scale atomistic simulations for such systems.


Introduction
A potential energy surface (PES) that yields potential energy of a system of atoms with given atomic coordinates is the fundamental enabler for atomistic simulation methods. In principle, ab initio or first principles methods that solve the Schrödinger equation, typically some approximation within the Kohn-Sham density functional theory (DFT) framework, 1,2 can be applied to directly calculate the PES. While such methods are highly accurate and transferable across diverse chemistries and bonding types, their high computational cost limit their application in molecular dynamics (MD) simulations to relatively small and simple systems containing up to a few hundreds of atoms and sub-nanosecond time scales. Empirical interatomic potentials, on the other hand, are a much cheaper alternative. The functional form of these potentials are drastically simplified with only a few fitting parameters to satisfy physical considerations. 3,4 However, the accuracy of the empirical potentials is necessarily limited by the approximations made in selecting the functional form, which are generally not transferable to another system with different bonding types.
In recent years, an alternative approach has gained popularity in constructing interatomic potentials with improved transferability. [5][6][7][8][9][10] In this approach, the atomic coordinates are featurized using local environment descriptors that are invariant to translations, rotations and permutations of homo-nuclear atoms, and are differentiable and unique. 8,11 A machine learning model is then trained to map the structural features to data (energies, forces, etc.) from first principles calculations. Such potentials have been demonstrated to achieve accuracy close to first principles methods at much lower computational costs. 5,7,9,10 The coefficients of the bispectrum of local atomic density were first applied in the Gaussian approximation potential by Bartók et al. 7 Thompson et al. later showed that a linear model of bispectrum coefficients from the lowest order -the so-called Spectral Neighbor Analysis Potential (SNAP) -can accurately reproduce DFT energies and forces as well as a variety of calculated properties (e.g., elastic constants and migration barrier for screw dislocations) in bcc Ta and W. 9,12 More recently, the current authors have extended the SNAP formalism to bcc Mo, fcc Ni and Cu, and the binary fcc Ni-bcc Mo alloy systems and showed that it outperforms traditional embedded atom method (EAM) and modified EAM potentials across a wide range of properties. 13,14 Thus far, SNAP models have mainly been developed for metallic systems. For ionic systems, a common strategy in constructing interatomic potentials is to incorporate long-ranged electrostatic interactions (e.g., through the use of the Ewald summation) on top of energy model. This has been done for both traditional empirical models 15,16 as well as modern local atomic environment descriptor-based potentials (e.g., GAP for the mixed ionic-covalent GaN 7 and neural network potential for ZnO 6 ). In this work, we develop a highly-accurate electrostatic SNAP (eSNAP) model for ionic α-Li 3 N (see Figure 1(a)). α-Li 3 N is one of the earliest lithium superionic conductors ever reported, 17

Optimized model parameters
In our proposed electrostatic SNAP (eSNAP) model, we write the total potential energy E p as the sum of the electrostatic contributions and the local (SNAP) energy due to the variations in atomic local environments (SNAP), as follows: where E el and E SN AP are the electrostatic energy computed using the Ewald summation approach 20 and the energy from SNAP, respectively, and γ is an effective screening prefactor for electrostatic interactions. An iterative procedure was developed to fit all model parameters using total energies and forces from DFT calculations until the training and test errors are converged (see Methods for details).
The final hyperparameters and coefficients for the optimized eSNAP model are given in Table 1. The optimized effective screening parameter γ is 0.057.    We also calculated the formation energy of Li Frenkel defects, which contributes significantly to Li + conductivity, in α-Li 3 N (Table 2). Kishida et al. 24 have shown earlier that vacancies are unlikely to form on Li1 sites, and our investigations have found that any Frenkel defect created by introducing a vacancy on Li1 always relaxes back to the pristine Li 3 N crystal in DFT. Here, we consider the two Frenkel configurations where a vacancy is introduced on a Li2 site and the interstitial Li is located at either the Li2 site (intraplanar, Figure 1

Structural Properties
or Li1 site (interplanar, Figure 1(c)). Both eSNAP and the Coulomb-Buckingham potential yield formation energies reasonably close to the DFT values, but both failed at reproducing the energy difference between the two configurations.
Finally, Figure 3 compares the calculated phonon dispersion curves of α-Li 3 N from eSNAP with those from DFT calculations. The phonon dispersion curves were calculated using the finite displacement approach on a 3 × 3 × 3 supercell as implemented in the phonopy package. 25 We find that the phonon dispersion curves calculated from eSNAP are in good agreement with that from DFT. The only discrepancy is the imaginary phonon mode at Γ point observed in DFT phonon dispersion. According to Wu et al. 26 , this lattice instability is associated with the vibration of Li2 sites along the c axis, resulting in a more stable phase that is only 0.3 meV/atom lower in energy after displacing Li2 site by ∼ 0.1Å. This energy difference is well within the energy prediction error of the eSNAP model. We also note that the experimentally measured phonon dispersion curves at room temperature do not exhibit this lattice instability. 23 In contrast, the phonon dispersion curve calculated from the Coulomb-Buckingham potential show severely overestimated frequencies ( Figure ??) due to its unsatisfactory force prediction.    These results indicate that grain boundaries may provide a rapid pathway for Li diffusion in α-Li 3 N. We also carried out MD simulation using the Coulomb-Buckingham potential, and the same trends are observed.

Discussion
In this work, we demonstrate that modern potentials based on local environment descriptors such as the Spectral Neighbor Analysis Potential (SNAP) can be adapted for ionic systems by incorporating long-range electrostatics. For highly ionic α-Li 3 N, we find that assigning formal charges of +1 and -3 to Li and N, respectively, with screening accounted for via a fitted hyperparameter (γ in Equation 1), result in a simpler, more stable potential model than variable charge models such as the charge equilibration (QEq) method. 29

Electrostatic SNAP (eSNAP) Model
The atomic environment around atom i at coordinates r can be described by its atomic neighbor density ρ i (r) with the following equation 7,9 : where r ii is the vector joining the coordinates of central atom i and its neighbor atom i , the cutoff function f c ensures that the neighbor atomic density decays smoothly to zero at cutoff radius R ii , and the dimensionless neighbor weights w i distinguish atoms of different types. This density function can be expanded as a generalized Fourier series in the 4D hyper-spherical harmonics U j m,m (θ, φ, θ 0 ) as follows: where the coefficients u j m,m are given by the inner product U j m,m |ρ . The bispectrum coefficients are then given as: where the constants H jmm j 1 m 1 m 1 j 2 m 2 m 2 are coupling coefficients.
In the original formulation of the non-ionic SNAP model, 9 the energy and forces are expressed as a linear function of the bispectrum coefficients, as follows: where α is the chemical identity of atoms and β α,k are the coefficients in the linear SNAP model for type α atoms.
R ii Interatomic distance local environment electrostatic interaction nuclei repulsion For ionic systems, electrostatic interactions spanning in the entire range of inter-atomic distances are indispensable in the construction of energy model due to the long-range tail beyond the cutoff distance for local environment description (see Figure 6). In our proposed electrostatic SNAP (eSNAP) model, we write the total potential energy as the sum of the electrostatic contributions and the local (SNAP) energy due to the variations in atomic local environments (SNAP), as follows: where E el is the electrostatic energy computed using the Ewald summation approach 20 and γ is an effective screening prefactor for electrostatic interactions. The coefficients (γ and β) can be solved by fitting the linear model to total energies and forces from DFT calculations.
In addition, nuclei repulsion emerge at extremely short inter-atomic distances. In this work, the Ziegler-Biersack-Littmark (ZBL) potential is used to account for short-ranged nuclei repulsion. 33 To ensure that the fitting process captures the relevant relationship between the bispectrum coefficients and the DFT energies and forces, the cutoff distances of ZBL were chosen to be short enough (R i = 1.0Å, R o = 1.5Å) such that the ZBL potential has negligible contribution to energies or forces among the initial training configurations where extremely close inter-atomic distances were not sampled.   Table 4 shows the weights applied on the different sets of training configurations during model training. As the initial training dataset contains many more configurations from AIMD snapshots with larger number of atoms, a much larger weight was applied on the energies of the distorted unit cells relative to those from the AIMD snapshots. A zero weight was applied on the negligibly small forces for the distorted unit cells. As shown in Figure 7a, the energies and forces differ greatly in magnitude and distribution due to differences in the scales and units. In the original SNAP training approach, the effect of this difference in magnitude and distribution is partially accounted for by treating the data weights as hyperparameters to be optimized. 13,14 In this work, we use the standardized z-scores of energies and forces (plotted in Figure 7b) as the targets in model training to avoid incorporating the effect of the distribution in the data weights, which are therefore fixed at the values in Table 4. The "standardized" eSNAP model in the fitting process is then given by the following:   e−ē σe . . .

Model Training and Test
where e is the energy per atom,ē is the mean of e, and σ e and σ F are the standard deviations of e and F , respectively. The mean of forces is omitted since it is close to zero. The coefficient vector β T to be solved can be written as: For bispectrum coefficient calculations, we used the implementation available in LAMMPS. 9 The two hyperparameters (cutoff distance R α and atomic weight w α ) for each element (Li and N in the case of Li 3 N) were determined using a two-step grid search scheme for the atomic weights and then followed by the cutoff distances. The mean absolute error (MAE) of forces from a linear model trained on the initial training set was chosen as the metric. For the atomic weights, it should be noted that the atomic density in ionic systems is generally higher than that in metallic systems; hence the search of atomic weights was performed in the range where |w α | < 1. Similarly, the search space for cutoff radius was limited to the range where R α < 4Å. The results from grid search (Figure ??) and the final hyperparameters (  procedure was repeated until there is no significant over-fitting. In this work, we use 150% of training MAE as the threshold to achieve a balance between the benefit gained by adding more training instances and the associated costs of performing more DFT calculations. It should be noted that this strategy is designed to bias the eSNAP model to improve the predictions on energy and force of MD snapshots, which is the target application of interest in this work.

DFT calculations
All DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) 35 within the projector augmented wave approach. 36 The Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation was adopted as the exchange-correlation functional. 37 . To ensure the convergence of energy and atomic force, a plane-wave energy cutoff of 520 eV and Γ-centered k-point meshes with a density of at least 30Å were employed for all static DFT calculations. For AIMD simulations, a single Γ k-point and a much lower energy cutoff of 300 eV were used for rapid propagation of trajectories.

Diffusivity calculations
The tracer diffusivity of Li D * is calculated from the mean square displacement (MSD) of all diffusing Li ions as described by the Einstein relation: where d is the number of dimensions in which diffusion occurs, N is the total number of diffusing Li ions, ∆r i (t) is the displacement of the ith Li ion at time t.
The charge diffusivity of Li D σ is calculated from the square net displacement of all diffusing Li ions, as described below: The Li conductivity at temperature T (unit: K) can be calculated from the charge diffusivity D σ using the Nernst-Einstein equation: where ρ is the molar density of Li, z is the charge of Li (+1), F is the Faraday constant, and R is the gas constant.
In addition, the ratio between the tracer and charge diffusivities is referred to in the literature as the Haven ratio H R = D * /D σ .
All the simulations with the eSNAP were performed using LAMMPS. 31 All the struc-