Study of Li atom diffusion in amorphous Li₃PO₄ with neural network potential Open Access

Atom/ion diffusion in amorphous materials is an important elementary process that has a great influence on the performance of various information and energy devices. So far, several theoretical studies have been carried out on the diffusion behaviors in various amorphous materials using first-principles methods.^1–6 However, the number of such studies is limited because of the following reasons. First, the diffusion pathway characterization in an amorphous structure is much more complicated than in crystals, due to the existence of many inequivalent pathways. Second, large supercells should be used to represent the stochastic nature of amorphous materials. Lastly, amorphous structures can often have variations in density and composition, the effects of which on the diffusion behaviors are desirable to be clarified.

Recently, a new method, which involves the construction of interatomic potentials via machine learning techniques to determine the relationship between atomic structures and corresponding energies, has attracted much attention as a promising way to achieve high computational accuracy and speed. In this method, interatomic potentials are “trained” using the results of first-principles calculations as the reference (or “training”) data for the learning process. Then it becomes capable of predicting structural energies that are not included in the training dataset. Such machine-learning-based interatomic potentials were found to achieve much higher accuracy than the conventional interatomic potentials with simple formulas. The machine-learning techniques used for this method include neural networks (NNs),^7–20 kernel regression,²¹ Gaussian process regression,²² linear regression,^23,24 and so on. Among these methods, NNs appear to be especially promising as demonstrated by Behler and co-workers.^25,26 They evaluated the performance of the NN potentials (NNPs) for various materials systems, such as Si,¹⁷ C,¹⁸ Cu,¹³ ZnO,¹⁵ TiO₂,²⁵ H₂O dimers,²⁷ Cu clusters supported on Zn oxide,¹⁴ and Au/Cu nanoparticles with water molecules.^26,28 Although the relatively large amount of required training data (up to ten thousands of reference DFT calculations) would be a weakness of NNs, their high accuracy is quite attractive.

Recently, the NN potential has been successfully applied to investigate the self-diffusion of surface defects on Cu low-index face.²⁹ However, to the best of our knowledge, so far there is no research where NN potentials were adopted to study atomic diffusion in amorphous materials. In order to extend the application of NN potentials to study this subject, a comprehensive validation of its reliability is necessary, considering that the atomic environment in amorphous materials is much more complicated than in crystals.

In the present research, we constructed NN potentials to study atom diffusion in amorphous materials and examined their reliability. The classical electrolyte, amorphous Li₃PO₄, was chosen as the benchmark material since it is one of the most widely used solid electrolytes for thin-film batteries and a number of other technologies.^30–38 The outline of this paper is as follows. Section II introduces the details on how to construct the NN potential for Li₃PO₄. The performance of the NN potential in energy prediction is examined in Sec. III. In Sec. IV, the performance of the NN potential on the Li diffusion behaviors in a disordered Li₃PO₄ model is discussed on the basis of our calculation using the NN potential together with the nudged elastic band, kinetic Monte Carlo, and molecular dynamics methods. The results are carefully compared with the DFT results. In Sec. V, using the NN potential, the construction of a large-scale model of amorphous Li₃PO₄, the structural features of the large model, and the Li diffusion behavior in it are examined. The conclusions are provided in Sec. VI.

All the reference DFT energies of training structures were calculated with Vienna ab initio simulation package (VASP).^39,40 The projector augmented wave (PAW) method was adopted to treat atomic core electrons, while the Perdew-Burke-Ernzerhof functional within the generalized gradient approximation was adopted to describe the electron-electron interactions.^41,42 The training structures have several sizes of supercells containing up to 64 atoms. The cutoff energy (for the plane wave basis set) and k-point mesh were set so that the convergence of the total energy within 1 meV/atom was achieved. Specifically, the cutoff energy was set to 700 eV, and the 4 × 4 × 4, 4 × 4 × 2, and 4 × 2 × 2 k-point meshes were adopted for the supercells having 15-16, 29-32 and 61-64 atoms, respectively.

The database for the training of the NN potential included the following three kinds of data: (1) Snapshots taken directly from the molecular dynamics (MD) trajectories at various temperatures from 300 K to 4000 K; (2) Defective structures created by just removing one Li atom or one Li₂O unit randomly from the MD snapshot structures; (3) Nudged-elastic-band (NEB) intermediate images obtained using a roughly trained NN potential. The structures of the third kind were included in the training dataset since we applied the NN potential to the investigation of diffusion paths. Since in an MD trajectory, a Li atom spends most of the time vibrating around equilibrium sites and only a small fraction of the time around the saddle points between the neighboring equilibrium sites, the first and second structures might contain very limited information related to transient structures around the saddle points which govern the diffusion behaviors. More specifically, the data of the third kind were obtained as follows: first, a rough NN potential was trained with the MD snapshots and defective structures; then, the NEB calculations were performed with the rough NN potential to obtain the Li vacancy diffusion paths in other amorphous structures. At last, the intermediate images, especially those close to the saddle points, were added into training database. The transient structures are of great significance in obtaining an accurate potential. The detailed discussion will be given in Sec. IV. The composition of the training database is given in the supplementary material (Table SIV).

In total, 38 592 structures were created. These structures were randomly distributed into the training and testing sets. The training set contains 30 874 (∼80%) structures, while the remaining 7718 (∼20%) structures were used as the independent testing data.

We adopted the NN potential proposed by Behler and his co-workers for its high accuracy and compatibility in multi-component systems.^{13–19,43–45} The basic idea of the NN potential is that the total energy of a structure can be expressed as the sum of energy contributions of respective atoms, and the each energy contribution can be predicted by the NN according to the local atomic environment of the atom. The following two kinds of atom-centered symmetry function (SF), the radial SF and angular SF, were used as the descriptors of the local atomic environment,

in which R_i is the interatomic distance between the central atom and atom i, θ_ij is the angle formed by the central atom at the vertex, atom i and atom j. f_c is the cutoff function that defines the local environment sphere,

where R_c is the cutoff distance, which is set to 7 Å in the present work. This value is determined on the basis of a simple convergence test, the results of which are shown in the supplementary material (Table SII). The properties of atom-centered SFs are elaborated in Refs. 46 and 47.

One may think that the long-range Coulomb interaction plays an important role in the Li₃PO₄, and its truncation at the cutoff radius might result in serious errors. In fact, Behler and his co-workers have proposed a high-dimensional NN scheme that can take account of the long-range Coulomb interaction explicitly.^15,27 In this method, a second set of atomic NNs has been constructed to predict the environment-dependent charges at respective atomic sites. They found that this method only marginally improves the accuracy in comparison with the original NN scheme, while the electrostatic energy part increases the complexity of the NN optimization and makes the fit more difficult.^45,48 Considering this, together with the fact that we succeeded in achieving reasonable accuracy for the Li₃PO₄ system without the above extension for the Coulomb interaction, we adopt the original Behler and Parrinello scheme of the NN potential in the present study.

In multi-component materials, the presence of different species should be treated separately. For instance, in Li₃PO₄, a Li atom can have other Li atoms, P atoms, and O atoms in its local atomic environment. Thus, for each species, we must employ three sets of radial SFs (for three possible neighboring elements) and six sets of angular SFs (note that for a central Li atom, the angular SF for the combination of O and P is equivalent to that of P and O). As can be seen from this, the number of SF types increases rapidly with the number of species. For this reason, the application of NN was restricted to systems with up to four elements so far.²⁸ 

Under the above situation, simplification of the atomic symmetry function set would be worth exploring. We tried this when we constructed the NN potential for Li₃PO₄. Our idea is to omit some types of angular SF that seem physically much less important. To be more specific, we considered the following features of bonding in the Li₃PO₄: P atoms tend to bond with four neighbor O atoms and form phosphates (PO₄³⁻) with the shape of tetrahedron, while Li and O atoms tend to form the Li₂O unit with covalent bonds. These P–O and Li–O covalent bonds have specific angular preferences. Based on these structure characteristics in Li₃PO₄, we used only four types of angular SF sets, i.e, P–O O, O–P P, Li–O O, and O–Li Li (the first letter denotes the element of the central vertex atom and the two letters after the dash denotes the elements of the neighbor atoms) while keeping all the radial SFs in constructing the NN potential. In fact, we examined several NN potentials constructed with different choices of angular SFs. In all the cases, the NN potentials were trained with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm,⁴² using the training set as reference data, while the testing set was used to validate them. Table I shows the lowest root mean squared errors (RMSEs) of energy prediction of the constructed NN potentials. Note that for each kind of NN potential, we repeated the training 10 times to get the best fitting. As seen in Table I, compared with the NN potential having no angular SF as inputs (No. 5), the NN potentials constructed using the angular SF sets have clearly improved accuracy. On the other hand, the difference in the degree of accuracy is hardly seen among NN potentials with different choices of angular SFs (No. 1 to No. 4). The results shown in Table I suggest that the description of the characteristics of P–O and Li–O bond angles is of great significance in constructing accurate NN potentials, while the description of the rest of angular SF sets is much less important. Meanwhile, using simplified input SF sets reduces the number of input nodes, which leads to the significant reduction of the total number of parameters in the NN potential.

Although this is just a preliminary attempt to simplify the NN potential, our idea successfully reduced the size of input SFs without severe loss in accuracy. We think this idea could be a promising solution to solve the dilemma of NN potential in the complex materials with many species. Of course, this approach has a limitation as well. Though in the case of Li₃PO₄, we can easily recognize the structural feature and corresponding important SF, this may not always be possible for other material systems. Another idea/approach should be explored in such cases if the simplification of NN potential is desirable.

On the basis of the comparison shown in Table I, we adopted the simplified input SF set No. 3 in the construction of the NN potential. The atomic NN for each species has two hidden layers with 15 hidden nodes in each layer. The hyperbolic tangent activation functions were employed at the nodes of the two hidden layers, and a linear activation function was employed at the output node. This NN architecture contains 3123 parameters in total.

In Fig. 1, the energies predicted by the NN potential were plotted against DFT energies. As mentioned above and shown in Table I, the residual mean square error (RMSE) of the energy prediction is 5.0 meV/atom for the training set and 5.6 meV/atom for the testing set. However, verifying the NN potential with the testing set might be problematic for the following two reasons. First, the structures from the testing set might be highly correlated with those in the training set due to the inclusion of continuous MD trajectories in the training database. Second, the prediction capacity of the NN in structures having much larger supercells is unclear since all the testing structures have relatively small supercells. Considering the above, we performed another verification test with data obtained using a totally different structure: 500 steps MD trajectory of Li₄₈P₁₆O₆₄ at 4000 K. The size of these testing structures is twice larger than the largest training structures, and these data are not correlated at all with any of the training structures. The comparison of DFT and NN potential energies was shown in Fig. 2. The RMSE of energy prediction is 6.2 meV/atom, which is similar to the error in the training and testing sets. This verification test further proved the accuracy of energy prediction and also demonstrated the transferability of the NN potential to larger structures.

First we applied the NN potential to the characterization of the two most important crystal phases of Li₃PO₄, which are labeled as β and γ. The crystalline γ-Li₃PO₄ has the orthorhombic Pnma structure, while the crystalline β-Li₃PO₄ has the orthorhombic Pmn2 structure. Table II lists the lattice constants obtained from DFT (GGA, PBE) calculations, NN potential calculations, and experiments. The lattice constants obtained by the NN potential are highly consistent with those by the DFT and experiments (the relative error does not exceed 1%). The total energies calculated by the NNP and DFT also agree with each other well. In agreement with DFT, the NN potential calculations show that the β form is energetically more favorable than the γ form. The energy difference per Li₃PO₄ unit between the the β and γ forms, E_γ − E_β, is 0.02 eV for the DFT and 0.01 eV for the NN potential.

In this section, we test the reliability of the NN potential in atomic diffusion behaviors. For this purpose, we examined diffusion behaviors by combining the NN potential with various computational approaches, i.e., NEB, kinetic Monte Carlo (KMC), and MD, and then compared the results with those obtained by the combination of these approaches with DFT. All computations were carried out using an amorphous Li₃PO₄ (a-Li₃PO₄) structure created with the typical melt quenching method: The γ-Li₃PO₄ structure (containing 32 atoms) was melted and equilibrated for 15 ps at 6000 K and then quenched to 300 K at the rate of 1 K/fs. At last, the structure was equilibrated for 3 ps at 300 K. The aforementioned process was performed with ab initio MD. This a-Li₃PO₄ structure was further relaxed with NN potential and DFT until the maximum force acting on atoms became less than 0.01 eV/Å to obtain the final models for the NN potential and DFT simulations, respectively.

The two final structures are almost the same: The discrepancy in the atomic positions between the two is 0.02 Å in average and 0.05 Å at a maximum. The total free energy is −200.93 eV in the model for the NN potential, while −200.90 eV for DFT. The final model for the DFT is shown in Fig. 3.

It should be noted that these models are partially crystallized after the annealing and relaxation procedures and thus may not be fully relevant to represent the structure of amorphous Li₃PO₄ synthesized experimentally via the pulsed laser deposition technique. Nevertheless, hereafter we call these models “amorphous” for simplicity. To construct a more reasonable model of amorphous Li₃PO₄ requires the use of much a larger supercell and a slightly different composition. The structure of amorphous Li₃PO₄ will be discussed in Sec. V.

The vacancy formation energy was defined as

in which E[V_Li] is the total energy of the supercell containing one Li vacancy, E[bulk] is the total energy of the disordered Li₃PO₄ supercell, and μ_Li is the chemical potential of Li obtained from DFT calculation of a perfect Li bcc crystal. The supercells with a vacancy have been generated by removing a Li atom and then relaxed the structure. Since the supercell contains 12 Li atoms that are not equivalent to one another, we can obtain 12 different Li vacancy formation energies. Figure 4 shows the comparison of the vacancy formation energies between DFT and NN potential. The average discrepancy is 0.029 eV and the maximum is 0.040 eV. The final vacancy structures relaxed by DFT and NN potential were also quite similar. The average and maximum differences in atomic positions are about 0.03 Å and 0.08 Å, respectively. From Fig. 4, we can also see that the NN potential systematically underestimated the energy of Li vacancy structures. Since the energies are predicted based on the structural similarity, it may occur that the NN potential gives a systematical error for one kind of configurations. However, due to the very complicated formula of NNs, it is difficult to analyze the origin of such an error. On the other hand, we would like to point out that the deviation is quite small (about 1 meV/atom) and less than the RMSE of energy prediction. So we do not think that this systematical underestimation affects the conclusion of the present study.

It should be noted that the large formation energies (∼5 eV) obtained here mean very small probability of formation of Li vacancy. In reality, a vacancy-interstitial pair would be formed rather than an isolated vacancy. Since the purpose of this section is the verification of the applicability of the NN potential in various calculations, we do not discuss the cases with complex defects here.

Here we examine the Li vacancy diffusion network (including the Li migration paths, barrier energies, and the topology connection of paths) using the NEB method. This method is frequently used to find migration paths and their energies. In this method, a number of intermediate images, which are constrained by imaginary “spring” force, are interposed between the starting and ending equilibrium positions of the target atom. By optimizing the atomic coordinates in each image according to the forces given by both imaginary spring and real potential energy surface, we can find the minimum energy path of atom migration between the two equilibrium sites. The optimization can be based on the potential energy surface provided by the DFT calculation or classical interatomic potential (e.g., NN potential).

In the present study, we assumed that a Li atom can directly hop to the neighbor vacancy site when the distance between the two sites is shorter than 3.5 Å. Then, we obtained a Li diffusion network comprising 46 individual migration paths in the disordered Li₃PO₄ model. Each of these paths was characterized by the NEB calculation (based on either NN potential or DFT). The location and shape of these migration paths are shown in Fig. 5. As we can see, the diffusion paths predicted by the two methods agree well with each other.

The barrier energy is defined as the energy difference between the activated and initial states. Since the atom can hop along the diffusion path in the forward and backward directions, each diffusion path has two diffusion barrier energies. Figure 6 gives the comparison of activation energies obtained by the NN potential and DFT. The average barrier energy is 0.58 eV in the DFT calculation while 0.57 eV in the NN potential. The prediction error of barrier energies is 0.048 eV in average, which is about 8% compared with the average barrier energy.

Here we would like to emphasize that to achieve a high accuracy in predicting barrier energies, it is necessary to include Li diffusion transient structures (introduced in Sec. II) into the training database. To prove it, another NN potential without using the transient structures is constructed for comparison. We raised the percentage of the training set to 92% so that the size of the training set for this potential (30 386) is similar to that of the original potential (30 874). We found that the accuracy of this potential in NEB calculation is worse than that of the original one (the average error of barrier energies is 0.073 eV), though the energy prediction accuracy is not bad (RMSE: 5.4 meV/atom for the training set and 5.6 meV/atom for the testing set). More detailed information is provided in the supplementary material.

To evaluate the realistic atomic diffusion coefficients as well as the effective activation energies for diffusion, KMC simulations are often performed using the information of the diffusion networks obtained from NEB calculations. So, we performed KMC simulations employing the two Li diffusion networks, obtained with DFT and NN potential respectively, to check the validity of the NN potential further. The probability for the Li atom hopping from the site i to the neighboring vacancy site j was determined using the harmonic transition state theory

where the prefactor k represents the jump frequency (typical the order of 10¹³ s⁻¹) and ΔE_ij is the barrier energy of Li vacancy migration computed with NEB. The outline of the KMC algorithm together with the values of parameters used in the present simulations is given in the supplementary material.

The KMC simulations were run at different temperatures from 500 K to 1500 K. For each simulation, 100 × 10⁶ MC events were performed. The diffusion coefficient D can be evaluated from the statistically averaged mean square displacement (MSD) of atoms on the basis of the Einstein relation

Figure 7 shows the temperature dependence of Li diffusion coefficients in the Li₃PO₄ model, and the overall diffusion obeys an ideal Arrhenius behavior. By fitting the plots with straight lines, the effective diffusion activation energy of Li vacancy diffusion is estimated as 0.602 eV for NN potential and 0.561 eV for DFT. The discrepancy between the two activation energies is 0.041 eV.

Molecular dynamics (MD) simulations are powerful to study atomic diffusion since they trace the time evolution of atomic motions, as well as the combination of NEB and KMC simulations. By analyzing the atomic trajectories obtained from constant-temperature MD simulations, we can compute the diffusion coefficient at the specific temperature.

We performed molecular dynamics with the Li₃PO₄ models (Li₁₁P₄O₁₆) with a Li vacancy at 800 K, 1000 K, 1200 K, and 1500 K. All these MD simulations used the same initial structure, which is an optimized amorphous vacancy model obtained in Sec. IV B. The simulation was performed with the canonical ensemble (NVT) using the Verlet algorithm, and the time step is 2 fs. For each MD run, the system was first equilibrated for 10 ps, and the data used to calculate the diffusion properties were accumulated for the following 50 ps. The averaged mean squared displacements (MSDs)⁴⁹ of Li atoms at different temperatures are shown in Fig. 8. According to the Einstein equation [Eq. (6)], the diffusion coefficient can be directly obtained from the slope of the MSD plot. From Fig. 8, we can see that MSDs of Li ions increase almost linearly with time at 1000 K, 1200 K, and 1500 K. On the other hand, the MSD plot of 800 K deviates from the linear relationship considerably. This abnormal behavior can be ascribed to the low Li diffusivity at low temperature. To be more specific, too few atomic hops occur during the simulation time, and the reliable macroscopic diffusion properties cannot be extracted from them. The longer MD simulation (1 ns, with NN potential), which is shown in the supplementary material, shows the normal linear MSD-time relationship. The Li diffusion coefficients extracted from the NN potential trajectories are 5.15 × 10⁻⁶ cm²/s at 1000 K, 1.17 × 10⁻⁵ cm²/s at 1200 K, and 3.20 × 10⁻⁵ cm²/s at 1500 K, while those from ab initio MD trajectories are 4.75 × 10⁻⁶ cm²/s at 1000 K, 1.01 × 10⁻⁵ cm²/s at 1200 K, and 3.17 × 10⁻⁵ cm²/s at 1500 K. The Li diffusion coefficients determined from the NN potential MD simulations are in good agreement with the DFT results.

In order to further confirm the agreement between the NN potential and DFT in the MD simulations, we calculated the total radial distribution functions averaged over the MD trajectories at respective temperatures, which are shown in Fig. 9. As can be seen in this figure, there are only slight discrepancies between the radial distribution functions obtained from the NN potential and DFT simulations.

The calculation speed of the NN potential MD is much faster than the ab initio MD. For example, generating a trajectory described above takes about 96 h on 192 cores (Intel^® Xeon^® processor L5640) in the case of ab initio method, while it takes only 4.5 h on a single core of the same processor in the case of the NN potential. That is, the calculation speed of NN potential MD simulations is about 4000 times faster than the ab initio ones.

Finally, we carried out long time MD simulations (1 ns per trajectory) to obtain reliable Li diffusivities. The results were plotted in Fig. 7 together with the KMC results, and the mean squared displacements of Li along these trajectories are shown in the supplementary material (Fig. S2). As can be seen from Fig. 7, the temperature dependence of the evaluated Li diffusivities obeys the Arrhenius law. The effective activation energy for diffusion obtained from the linear fitting of the Arrhenius law is 0.59 eV for the MD results with the NN potential, which qualitatively agrees with our KMC simulation results using the same NN potential. The quantitative difference in the diffusion coefficients seen in Fig. 7 between the MD and KMC can be attributed to the jump frequency used in our KMC simulations: The same value of 10¹³ s⁻¹ was set to the jump frequencies of all the transitions, which must be a very crude approximation.

By means of DFT calculation, the structure and properties of crystalline Li₃PO₄ and LiPON have been studied,^32,34,50,51 but the theoretical research on amorphous Li₃PO₄ is still scarce. We used the NN potential developed in the present work to construct large-scale amorphous Li₃PO₄ models. According to the experimental observation,^38,52 the composition of a Li₃PO₄ thin film fabricated by pulsed laser deposition is often slightly different from the stoichiometric one. The ratio of Li/P is about 2.9 according to the inductively coupled plasma atomic emission spectroscopy.³⁸ The deviation from the stoichiometric value can be attributed to the partial condensation reaction of Li₃PO₄, i.e., 2Li₃PO₄ → Li₄P₂O₇ + Li₂O, where the resultant Li₂O is lost during the fabrication.

The amorphous models having similar Li/P ratio were generated as follows. We started from three crystalline supercells composed of 16, 64, and 128 Li₃PO₄ units. Then, one, three, and six Li₂O units were removed from the supercells to set the Li/P ratio as 2.875, 2.906, and 2.906, respectively. The supercells contain 125, 503, and 1006 atoms, respectively. For the three initial structures, the amorphous models were generated by the simulated annealing method with the NN-potential-based MD. For comparison, ab initio MD was also used to perform the same simulated annealing procedures on the smallest model (Li₄₆P₁₆O₆₃). The detailed procedure is as follows: (1) initial structures were heated to 2000 K and thermalized for 30 ps; (2) the structures were subsequently quenched to 300 K with a speed of 0.5 K/fs; (3) the final structures were equilibrated for 5 ps at 300 K and then relaxed until the maximum force acting on an atom was smaller than 0.01 eV/Å.

The resultant structures are shown in Fig. 10. We can see a strong similarity between the structures of Li₄₆P₁₆O₆₃ created by the NN potential and ab initio MD. This demonstrates that the NN potential can reproduce final structures of the melt-quench DFT simulations. On the other hand, it is also noticeable that the generated Li₄₆P₁₆O₆₃ structures are partially crystallized. The arrangements of Li and P atoms in the Li₄₆P₁₆O₆₃ structures are analogous to those in γ and β-Li₃PO₄, though the orientations of P–O tetrahedrons are disordered. This partial crystallization problem was eliminated by using the larger supercells. As seen in Fig. 10(a), the resultant structures of large models (Li₁₈₆P₆₄O₂₅₃ and Li₃₇₂P₁₂₈O₅₀₆) are much more disordered. The broader peaks in radial distribution functions shown in Fig. 10(b) also suggest this. It should be stressed that all the resultant structures include not only the isolated PO₄³⁻ tetrahedrons but also the P₂O₇⁴⁻ dimers. The P₂O₇⁴⁻ dimers are formed from two PO₄ tetrahedrons via corner-sharing. The shape of P₂O₇⁴⁻ dimers found inside the resultant structures is shown in Fig. 10(a). The presence of this dimer structure is consistent with the observed Fourier transform infrared spectroscopy of Li₃PO₄ thin film.³⁸ 

The Li diffusion coefficients in the largest amorphous Li₃PO₄ model were evaluated by MD simulations for 100 ps at 600 K, 700 K, 800 K, 1000 K, and 1200 K. The obtained Li diffusion coefficients are plotted in Fig. 11. The activation energy obtained from the fitting to the Arrhenius law is 0.55 eV.

We also plotted Li diffusion coefficients that were measured experimentally by Kuwata et al. in Fig. 11. The Li diffusion coefficient D_Li was measured by SIMS for ⁶Li/⁷Li diffusion couples prepared by an ion-exchange method using liquid electrolyte and a mask method based on thin-film deposition. The ionic conductivity σ was measured with impedance spectroscopy, and then the diffusion coefficient D_σ is calculated with the Nernst-Einstein equation. The diffusion activation energies obtained from these experiments are 0.58 eV (SIMS, ion-exchange), 0.57 eV (SIMS, mask), and 0.55 eV (impedance spectroscopy). These values agree well with the value estimated based on our simulation, 0.55 eV.

We investigated the application of the high-dimensional neural network (NN) potential in the study of atomic diffusion using Li₃PO₄ as the model system. For this purpose, an NN potential was constructed using about 38 592 various structures and corresponding energies. The constructed NN potential is capable of accurately predicting the energy of these structures, and the structural and energetic properties of crystalline Li₃PO₄ calculated with the NN potential are in good agreement with DFT calculation and experiment. In addition, the effective method to reduce the number of input symmetry function descriptors was discussed. Such a method might be useful in applying the NN potential to systems with more chemical elements than Li₃PO₄.

The application of NN potentials in atomic diffusion was explored by computing Li vacancy formation energies, diffusion pathways and barrier energies, and diffusion coefficients and activation energies in a small disordered Li₃PO₄ model. We found that the NN potential can accurately reproduce the DFT results in all of the above calculations with much smaller computational costs. Further, it was shown that a large-scale amorphous Li₃PO₄ structure can be constructed by the NN potential. The Li diffusivity in it was determined by MD simulation and agrees well with the experimental observation.

See supplementary material for the parameters of symmetry functions, the details of kinetic Monte Carlo and molecular dynamics simulations, and the performance of potential trained without transient structures.

This work is partially supported by JSPS KAKENHI Grant Number JP16H01535, Japan, the Support Program for Starting up Innovation hub from Japan Science and Technology Agency (JST), and CREST-JST, Japan. Part of the computations in this work was performed using the facilities of the Supercomputer Center, the Institute for Solid State Physics, and the University of Tokyo.