To clarify atomic diffusion in amorphous materials, which is important in novel information and energy devices, theoretical methods having both reliability and computational speed are eagerly anticipated. In the present study, we applied neural network (NN) potentials, a recently developed machine learning technique, to the study of atom diffusion in amorphous materials, using Li_{3}PO_{4} as a benchmark material. The NN potential was used together with the nudged elastic band, kinetic Monte Carlo, and molecular dynamics methods to characterize Li vacancy diffusion behavior in the amorphous Li_{3}PO_{4} model. By comparing these results with corresponding DFT calculations, we found that the average error of the NN potential is 0.048 eV in calculating energy barriers of diffusion paths, and 0.041 eV in diffusion activation energy. Moreover, the diffusion coefficients obtained from molecular dynamics are always consistent with those from *ab initio* molecular dynamics simulation, while the computation speed of the NN potential is 3–4 orders of magnitude faster than DFT. Lastly, the structure of amorphous Li_{3}PO_{4} and the ion transport properties in it were studied with the NN potential using a large supercell model containing more than 1000 atoms. The formation of P_{2}O_{7} units was observed, which is consistent with the experimental characterization. The Li diffusion activation energy was estimated to be 0.55 eV, which agrees well with the experimental measurements.

## I. INTRODUCTION

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,^{21} Gaussian process regression,^{22} 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,^{17} C,^{18} Cu,^{13} ZnO,^{15} TiO_{2},^{25} H_{2}O dimers,^{27} Cu clusters supported on Zn oxide,^{14} 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.^{29} 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_{3}PO_{4}, 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_{3}PO_{4}. 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_{3}PO_{4} 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_{3}PO_{4}, the structural features of the large model, and the Li diffusion behavior in it are examined. The conclusions are provided in Sec. VI.

## II. METHOD

### A. DFT calculations for reference energies

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_{2}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.

### B. Construction of neural network potential

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_{3}PO_{4}, 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_{3}PO_{4} 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.

### C. Simplification of neural network potential

In multi-component materials, the presence of different species should be treated separately. For instance, in Li_{3}PO_{4}, 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.^{28}

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_{3}PO_{4}. 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_{3}PO_{4}: P atoms tend to bond with four neighbor O atoms and form phosphates (PO_{4}^{3−}) with the shape of tetrahedron, while Li and O atoms tend to form the Li_{2}O unit with covalent bonds. These P–O and Li–O covalent bonds have specific angular preferences. Based on these structure characteristics in Li_{3}PO_{4}, 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,^{42} 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.

. | . | . | . | . | RMSE (meV/atom) . | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|

. | Species . | Types of angular symmetry function set . | Network . | No. of parameters . | Training . | Testing . | |||||

Li | Li–LiLi | Li–PP | Li–OO | Li–LiP | Li–LiO | Li–PO | 156-15-15-1 | 7833 | 4.9 | 5.7 | |

1 | P | P–LiLi | P–PP | P–OO | P–LiP | P–LiO | P–PO | 156-15-15-1 | |||

O | O–LiLi | O–PP | O–OO | O–LiP | O–LiO | O–PO | 156-15-15-1 | ||||

Li | Li–LiLi | Li–PP | Li–OO | 90-15-15-1 | 4863 | 5.0 | 5.5 | ||||

2 | P | P–LiLi | P–PP | P–OO | 90-15-15-1 | ||||||

O | O–LiLi | O–PP | O–OO | 90-15-15-1 | |||||||

Li | Li–OO | 46-15-15-1 | 3213 | 5.0 | 5.6 | ||||||

3 | P | P–OO | 46-15-15-1 | ||||||||

O | O–LiLi | O–PP | 68-15-15-1 | ||||||||

Li | 24-15-15-1 | 2553 | 5.7 | 6.0 | |||||||

4 | P | P–OO | 46-15-15-1 | ||||||||

O | O–PP | 46-15-15-1 | |||||||||

Li | 24-15-15-1 | 1893 | 8.3 | 8.6 | |||||||

5 | P | 24-15-15-1 | |||||||||

O | 24-15-15-1 |

. | . | . | . | . | RMSE (meV/atom) . | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|

. | Species . | Types of angular symmetry function set . | Network . | No. of parameters . | Training . | Testing . | |||||

Li | Li–LiLi | Li–PP | Li–OO | Li–LiP | Li–LiO | Li–PO | 156-15-15-1 | 7833 | 4.9 | 5.7 | |

1 | P | P–LiLi | P–PP | P–OO | P–LiP | P–LiO | P–PO | 156-15-15-1 | |||

O | O–LiLi | O–PP | O–OO | O–LiP | O–LiO | O–PO | 156-15-15-1 | ||||

Li | Li–LiLi | Li–PP | Li–OO | 90-15-15-1 | 4863 | 5.0 | 5.5 | ||||

2 | P | P–LiLi | P–PP | P–OO | 90-15-15-1 | ||||||

O | O–LiLi | O–PP | O–OO | 90-15-15-1 | |||||||

Li | Li–OO | 46-15-15-1 | 3213 | 5.0 | 5.6 | ||||||

3 | P | P–OO | 46-15-15-1 | ||||||||

O | O–LiLi | O–PP | 68-15-15-1 | ||||||||

Li | 24-15-15-1 | 2553 | 5.7 | 6.0 | |||||||

4 | P | P–OO | 46-15-15-1 | ||||||||

O | O–PP | 46-15-15-1 | |||||||||

Li | 24-15-15-1 | 1893 | 8.3 | 8.6 | |||||||

5 | P | 24-15-15-1 | |||||||||

O | 24-15-15-1 |

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_{3}PO_{4}, 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.

## III. SIMPLE TEST FOR NEURAL NETWORK POTENTIAL

### A. Energy prediction

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_{48}P_{16}O_{64} 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.

### B. Crystal structure parameters

First we applied the NN potential to the characterization of the two most important crystal phases of Li_{3}PO_{4}, which are labeled as *β* and *γ*. The crystalline *γ*-Li_{3}PO_{4} has the orthorhombic *Pnma* structure, while the crystalline *β*-Li_{3}PO_{4} 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_{3}PO_{4} unit between the the *β* and *γ* forms, *E*_{γ} − *E*_{β}, is 0.02 eV for the DFT and 0.01 eV for the NN potential.

. | NNP . | DFT . | Expt.^{53,54}
. | ||
---|---|---|---|---|---|

. | Energy (eV) . | Lattice constant (Å) . | Energy (eV) . | Lattice constant (Å) . | Lattice constant (Å) . |

γ- Li_{3}PO_{4} | −200.963 | a = 5.003 | −200.935 | a = 5.001 | a = 4.927 |

b = 6.108 | b = 6.177 | b = 6.120 | |||

c = 10.502 | c = 10.615 | c = 10.490 | |||

β- Li_{3}PO_{4} | −100.491 | a = 4.914 | −100.484 | a = 4.923 | a = 4.856 |

b = 5.287 | b = 5.298 | b = 5.240 | |||

c = 6.164 | c = 6.177 | c = 6.115 |

. | NNP . | DFT . | Expt.^{53,54}
. | ||
---|---|---|---|---|---|

. | Energy (eV) . | Lattice constant (Å) . | Energy (eV) . | Lattice constant (Å) . | Lattice constant (Å) . |

γ- Li_{3}PO_{4} | −200.963 | a = 5.003 | −200.935 | a = 5.001 | a = 4.927 |

b = 6.108 | b = 6.177 | b = 6.120 | |||

c = 10.502 | c = 10.615 | c = 10.490 | |||

β- Li_{3}PO_{4} | −100.491 | a = 4.914 | −100.484 | a = 4.923 | a = 4.856 |

b = 5.287 | b = 5.298 | b = 5.240 | |||

c = 6.164 | c = 6.177 | c = 6.115 |

## IV. PERFORMANCE OF NEURAL NETWORK POTENTIAL IN STUDYING ATOMIC DIFFUSION

### A. Structural model for this study

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_{3}PO_{4} (*a*-Li_{3}PO_{4}) structure created with the typical melt quenching method: The *γ*-Li_{3}PO_{4} 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_{3}PO_{4} 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_{3}PO_{4} 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_{3}PO_{4} requires the use of much a larger supercell and a slightly different composition. The structure of amorphous Li_{3}PO_{4} will be discussed in Sec. V.

### B. Vacancy formation energy

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_{3}PO_{4} 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.

### C. Diffusion paths and network

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_{3}PO_{4} 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.

### D. Kinetic Monte Carlo simulation

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^{13} s^{−1}) 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^{6} 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_{3}PO_{4} 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.

### E. Molecular dynamics simulation

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_{3}PO_{4} models (Li_{11}P_{4}O_{16}) 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)^{49} 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^{−6} cm^{2}/s at 1000 K, 1.17 × 10^{−5} cm^{2}/s at 1200 K, and 3.20 × 10^{−5} cm^{2}/s at 1500 K, while those from *ab initio* MD trajectories are 4.75 × 10^{−6} cm^{2}/s at 1000 K, 1.01 × 10^{−5} cm^{2}/s at 1200 K, and 3.17 × 10^{−5} cm^{2}/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^{13} s^{−1} was set to the jump frequencies of all the transitions, which must be a very crude approximation.

## V. LARGE-SCALE SIMULATION OF Li DIFFUSION IN AMORPHOUS Li_{3}PO_{4}

### A. Amorphous Li_{3}PO_{4} structure

By means of DFT calculation, the structure and properties of crystalline Li_{3}PO_{4} and LiPON have been studied,^{32,34,50,51} but the theoretical research on amorphous Li_{3}PO_{4} is still scarce. We used the NN potential developed in the present work to construct large-scale amorphous Li_{3}PO_{4} models. According to the experimental observation,^{38,52} the composition of a Li_{3}PO_{4} 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.^{38} The deviation from the stoichiometric value can be attributed to the partial condensation reaction of Li_{3}PO_{4}, i.e., 2Li_{3}PO_{4} → Li_{4}P_{2}O_{7} + Li_{2}O, where the resultant Li_{2}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_{3}PO_{4} units. Then, one, three, and six Li_{2}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_{46}P_{16}O_{63}). 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_{46}P_{16}O_{63} 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_{46}P_{16}O_{63} structures are partially crystallized. The arrangements of Li and P atoms in the Li_{46}P_{16}O_{63} structures are analogous to those in γ and β-Li_{3}PO_{4}, 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_{186}P_{64}O_{253} and Li_{372}P_{128}O_{506}) 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_{4}^{3−} tetrahedrons but also the P_{2}O_{7}^{4−} dimers. The P_{2}O_{7}^{4−} dimers are formed from two PO_{4} tetrahedrons via corner-sharing. The shape of P_{2}O_{7}^{4−} 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_{3}PO_{4} thin film.^{38}

### B. Diffusivity of Li

The Li diffusion coefficients in the largest amorphous Li_{3}PO_{4} 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 ^{6}Li/^{7}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.

## VI. CONCLUSIONS

We investigated the application of the high-dimensional neural network (NN) potential in the study of atomic diffusion using Li_{3}PO_{4} 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_{3}PO_{4} 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_{3}PO_{4}.

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_{3}PO_{4} 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_{3}PO_{4} 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.