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 Li3PO4 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 Li3PO4 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 Li3PO4 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 P2O7 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.

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 TiO2,25 H2O 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 Li3PO4, 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 Li3PO4. 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 Li3PO4 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 Li3PO4, 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 Li2O 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,

G1=ieηRi2fc(Ri),
(1)
G2=21ζiji1+λcosθij×eηRi2+Rj2fc(Ri)fc(Rj),
(2)

in which Ri 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. fc is the cutoff function that defines the local environment sphere,

fcRij=0.5×cosπRiRc+1Ri<Rc0Ri>Rc,
(3)

where Rc 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 Li3PO4, 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 Li3PO4 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 Li3PO4, 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 Li3PO4. 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 Li3PO4: P atoms tend to bond with four neighbor O atoms and form phosphates (PO43−) with the shape of tetrahedron, while Li and O atoms tend to form the Li2O unit with covalent bonds. These P–O and Li–O covalent bonds have specific angular preferences. Based on these structure characteristics in Li3PO4, 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.

TABLE I.

The root mean squared errors (RMSEs) of the NN energies obtained for the training set and the testing set of the Li3PO4 system using different input symmetry function sets. The number of parameters in the NN potential is also given for each architecture.

RMSE (meV/atom)
SpeciesTypes of angular symmetry function setNetworkNo. of parametersTrainingTesting
 Li Li–LiLi Li–PP Li–OO Li–LiP Li–LiO Li–PO 156-15-15-1 7833 4.9 5.7 
P–LiLi P–PP P–OO P–LiP P–LiO P–PO 156-15-15-1 
 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 
P–LiLi P–PP P–OO    90-15-15-1 
 O–LiLi O–PP O–OO    90-15-15-1 
 Li Li–OO      46-15-15-1 3213 5.0 5.6 
P–OO      46-15-15-1 
 O–LiLi O–PP     68-15-15-1 
 Li       24-15-15-1 2553 5.7 6.0 
P–OO      46-15-15-1 
 O–PP      46-15-15-1 
 Li       24-15-15-1 1893 8.3 8.6 
      24-15-15-1 
       24-15-15-1 
RMSE (meV/atom)
SpeciesTypes of angular symmetry function setNetworkNo. of parametersTrainingTesting
 Li Li–LiLi Li–PP Li–OO Li–LiP Li–LiO Li–PO 156-15-15-1 7833 4.9 5.7 
P–LiLi P–PP P–OO P–LiP P–LiO P–PO 156-15-15-1 
 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 
P–LiLi P–PP P–OO    90-15-15-1 
 O–LiLi O–PP O–OO    90-15-15-1 
 Li Li–OO      46-15-15-1 3213 5.0 5.6 
P–OO      46-15-15-1 
 O–LiLi O–PP     68-15-15-1 
 Li       24-15-15-1 2553 5.7 6.0 
P–OO      46-15-15-1 
 O–PP      46-15-15-1 
 Li       24-15-15-1 1893 8.3 8.6 
      24-15-15-1 
       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 Li3PO4, 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 Li48P16O64 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.

FIG. 1.

Comparison between the density functional theory (DFT) and NN potential on the total energies of the (a) training set and (b) testing set structures.

FIG. 1.

Comparison between the density functional theory (DFT) and NN potential on the total energies of the (a) training set and (b) testing set structures.

Close modal
FIG. 2.

Comparison of the total energy changes along molecular dynamics trajectories between the DFT and NN potential (NNP). The models used in this simulation contain 128 atoms, and none of the structures in the trajectories were included in the training set. The residual mean square error (RMSE) of NNP is 6.2 meV/atom.

FIG. 2.

Comparison of the total energy changes along molecular dynamics trajectories between the DFT and NN potential (NNP). The models used in this simulation contain 128 atoms, and none of the structures in the trajectories were included in the training set. The residual mean square error (RMSE) of NNP is 6.2 meV/atom.

Close modal

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

TABLE II.

Comparison of the lattice parameters and total energies between β- and γ-Li3PO4 crystal structures obtained with NN potential, DFT, and experimental measurement.

NNPDFTExpt.53,54
Energy (eV)Lattice constant (Å)Energy (eV)Lattice constant (Å)Lattice constant (Å)
γ- Li3PO4 −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 
β- Li3PO4 −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 
NNPDFTExpt.53,54
Energy (eV)Lattice constant (Å)Energy (eV)Lattice constant (Å)Lattice constant (Å)
γ- Li3PO4 −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 
β- Li3PO4 −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 

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 Li3PO4 (a-Li3PO4) structure created with the typical melt quenching method: The γ-Li3PO4 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-Li3PO4 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.

FIG. 3.

The final structure of the amorphous Li3PO4 model (Li12P4O16) after melt-quenching simulation. The structure was fully optimized with DFT.

FIG. 3.

The final structure of the amorphous Li3PO4 model (Li12P4O16) after melt-quenching simulation. The structure was fully optimized with DFT.

Close modal

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 Li3PO4 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 Li3PO4 requires the use of much a larger supercell and a slightly different composition. The structure of amorphous Li3PO4 will be discussed in Sec. V.

The vacancy formation energy was defined as

Ef=E[VLi]E[bulk]+μLi,
(4)

in which E[VLi] is the total energy of the supercell containing one Li vacancy, E[bulk] is the total energy of the disordered Li3PO4 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.

FIG. 4.

Li vacancy formation energies in the amorphous Li12P4O16 model calculated by the DFT and NN potential.

FIG. 4.

Li vacancy formation energies in the amorphous Li12P4O16 model calculated by the DFT and NN potential.

Close modal

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 Li3PO4 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.

FIG. 5.

All diffusion paths of Li atoms in the amorphous Li12P4O16 model, predicted by DFT and NN potential with the nudged elastic band method.

FIG. 5.

All diffusion paths of Li atoms in the amorphous Li12P4O16 model, predicted by DFT and NN potential with the nudged elastic band method.

Close modal

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.

FIG. 6.

Comparison of barrier energies of 92 vacancy diffusion paths in the Li3PO4 model calculated by the DFT and NN potential. (a) The barrier energies obtained by the DFT and NN potential and (b) the difference between them.

FIG. 6.

Comparison of barrier energies of 92 vacancy diffusion paths in the Li3PO4 model calculated by the DFT and NN potential. (a) The barrier energies obtained by the DFT and NN potential and (b) the difference between them.

Close modal

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

p=kexpΔEijkBT,
(5)

where the prefactor k represents the jump frequency (typical the order of 1013 s−1) and ΔEij 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 × 106 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

D=16limΔtMSD=16limΔtrit0+Δtrit02Δt.
(6)

Figure 7 shows the temperature dependence of Li diffusion coefficients in the Li3PO4 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.

FIG. 7.

Temperature dependence of Li diffusion coefficients obtained from the kinetic Monte Carlo and molecular dynamics simulations. The diffusion network topology and hopping possibilities obtained from the NEB calculations (for details, see Sec. III C) were used as the parameters for the KMC simulation.

FIG. 7.

Temperature dependence of Li diffusion coefficients obtained from the kinetic Monte Carlo and molecular dynamics simulations. The diffusion network topology and hopping possibilities obtained from the NEB calculations (for details, see Sec. III C) were used as the parameters for the KMC simulation.

Close modal

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 Li3PO4 models (Li11P4O16) 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 cm2/s at 1000 K, 1.17 × 10−5 cm2/s at 1200 K, and 3.20 × 10−5 cm2/s at 1500 K, while those from ab initio MD trajectories are 4.75 × 10−6 cm2/s at 1000 K, 1.01 × 10−5 cm2/s at 1200 K, and 3.17 × 10−5 cm2/s at 1500 K. The Li diffusion coefficients determined from the NN potential MD simulations are in good agreement with the DFT results.

FIG. 8.

Average mean square displacements along the molecular dynamics trajectories at 800 K, 1000 K, 1200 K, and 1500 K. The MSDs are averaged for 15 ps over 50 ps trajectories.

FIG. 8.

Average mean square displacements along the molecular dynamics trajectories at 800 K, 1000 K, 1200 K, and 1500 K. The MSDs are averaged for 15 ps over 50 ps trajectories.

Close modal

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.

FIG. 9.

Radial distribution functions of Li3PO4 averaged over MD trajectories at 800, 100, 1200, and 1500 K.

FIG. 9.

Radial distribution functions of Li3PO4 averaged over MD trajectories at 800, 100, 1200, and 1500 K.

Close modal

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 1013 s−1 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 Li3PO4 and LiPON have been studied,32,34,50,51 but the theoretical research on amorphous Li3PO4 is still scarce. We used the NN potential developed in the present work to construct large-scale amorphous Li3PO4 models. According to the experimental observation,38,52 the composition of a Li3PO4 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 Li3PO4, i.e., 2Li3PO4 → Li4P2O7 + Li2O, where the resultant Li2O 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 Li3PO4 units. Then, one, three, and six Li2O 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 (Li46P16O63). 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 Li46P16O63 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 Li46P16O63 structures are partially crystallized. The arrangements of Li and P atoms in the Li46P16O63 structures are analogous to those in γ and β-Li3PO4, 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 (Li186P64O253 and Li372P128O506) 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 PO43− tetrahedrons but also the P2O74− dimers. The P2O74− dimers are formed from two PO4 tetrahedrons via corner-sharing. The shape of P2O74− 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 Li3PO4 thin film.38 

FIG. 10.

The structures of amorphous Li3PO4 created by the melt-quenching simulations. (a) The resultant structures with different size, and the shape of P2O74− dimmers inside them. (b) The radial distribution functions of corresponding amorphous structures. Each function is averaged over 5 ps MD run at 300 K.

FIG. 10.

The structures of amorphous Li3PO4 created by the melt-quenching simulations. (a) The resultant structures with different size, and the shape of P2O74− dimmers inside them. (b) The radial distribution functions of corresponding amorphous structures. Each function is averaged over 5 ps MD run at 300 K.

Close modal

The Li diffusion coefficients in the largest amorphous Li3PO4 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.

FIG. 11.

The diffusivity of Li in amorphous Li3PO4. The simulation results were obtained from the large-scale MD (with the supercell of 1006 atoms) using NN potentials. The experimental results were measured with different methods (taken from Ref. 52).

FIG. 11.

The diffusivity of Li in amorphous Li3PO4. The simulation results were obtained from the large-scale MD (with the supercell of 1006 atoms) using NN potentials. The experimental results were measured with different methods (taken from Ref. 52).

Close modal

We also plotted Li diffusion coefficients that were measured experimentally by Kuwata et al. in Fig. 11. The Li diffusion coefficient DLi was measured by SIMS for 6Li/7Li 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 Li3PO4 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 Li3PO4 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 Li3PO4.

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 Li3PO4 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 Li3PO4 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.

1.
S.
Hao
and
C.
Wolverton
,
J. Phys. Chem. C
117
,
8009
(
2013
).
2.
P.
Johari
,
Y.
Qi
, and
V. B.
Shenoy
,
Nano Lett.
11
,
5494
(
2011
).
3.
D. M.
Guzman
,
N.
Onofrio
, and
A.
Strachan
,
J. Appl. Phys.
117
,
195702
(
2015
).
4.
G. A.
Tritsaris
,
K.
Zhao
,
O. U.
Okeke
, and
E.
Kaxiras
,
J. Phys. Chem. C
116
,
22212
(
2012
).
5.
N.
Onofrio
,
D.
Guzman
, and
A.
Strachan
,
Nat. Mater.
14
,
440
(
2015
).
6.
K.
Sankaran
,
L.
Goux
,
S.
Clima
,
M.
Mees
,
J. A.
Kittl
,
M.
Jurczak
,
L.
Altimime
,
G.-M.
Rignanese
, and
G.
Pourtois
,
ECS Trans.
45
,
317
(
2012
).
7.
K. V. J.
Jose
,
N.
Artrith
, and
J.
Behler
,
J. Chem. Phys.
136
,
194111
(
2012
).
8.
M.
Malshe
,
R.
Narulkar
,
L. M.
Raff
,
M.
Hagan
,
S.
Bukkapatnam
, and
R.
Komanduri
,
J. Chem. Phys.
129
,
44111
(
2008
).
9.
M.
Malshe
,
R.
Narulkar
,
L. M.
Raff
,
M.
Hagan
,
S.
Bukkapatnam
,
P. M.
Agrawal
, and
R.
Komanduri
,
J. Chem. Phys.
130
,
184102
(
2009
).
10.
J.
Behler
,
S.
Lorenz
, and
K.
Reuter
,
J. Chem. Phys.
127
,
14705
(
2007
).
11.
J.
Ludwig
and
D. G.
Vlachos
,
J. Chem. Phys.
127
,
154716
(
2007
).
12.
A.
Pukrittayakamee
,
M.
Malshe
,
M.
Hagan
,
L. M.
Raff
,
R.
Narulkar
,
S.
Bukkapatnum
, and
R.
Komanduri
,
J. Chem. Phys.
130
,
134101
(
2009
).
13.
N.
Artrith
and
J.
Behler
,
Phys. Rev. B
85
,
045439
(
2012
).
14.
N.
Artrith
,
B.
Hiller
, and
J.
Behler
,
Phys. Status Solidi
250
,
1191
(
2013
).
15.
N.
Artrith
,
T.
Morawietz
, and
J.
Behler
,
Phys. Rev. B
83
,
153101
(
2011
).
16.
J.
Behler
,
R.
Martoňák
,
D.
Donadio
, and
M.
Parrinello
,
Phys. Rev. Lett.
100
,
185501
(
2008
).
17.
J.
Behler
and
M.
Parrinello
,
Phys. Rev. Lett.
98
,
146401
(
2007
).
18.
R. Z.
Khaliullin
,
H.
Eshet
,
T. D.
Kühne
,
J.
Behler
, and
M.
Parrinello
,
Nat. Mater.
10
,
693
(
2011
).
19.
R. Z.
Khaliullin
,
H.
Eshet
,
T. D.
Kühne
,
J.
Behler
, and
M.
Parrinello
,
Phys. Rev. B
81
,
100103
(
2010
).
20.
J.
Behler
,
J. Chem. Phys.
145
,
170901
(
2016
).
21.
M.
Rupp
,
A.
Tkatchenko
,
K.-R.
Müller
, and
O. A.
von Lilienfeld
,
Phys. Rev. Lett.
108
,
058301
(
2012
).
22.
A. P.
Bartók
,
M. C.
Payne
,
R.
Kondor
, and
G.
Csányi
,
Phys. Rev. Lett.
104
,
136403
(
2010
).
23.
A.
Seko
,
A.
Takahashi
, and
I.
Tanaka
,
Phys. Rev. B
90
,
024101
(
2014
).
24.
A.
Seko
,
A.
Takahashi
, and
I.
Tanaka
,
Phys. Rev. B
92
,
54113
(
2015
).
25.
N.
Artrith
and
A.
Urban
,
Comput. Mater. Sci.
114
,
135
(
2016
).
26.
N.
Artrith
and
A. M.
Kolpak
,
Comput. Mater. Sci.
110
,
20
(
2015
).
27.
T.
Morawietz
,
V.
Sharma
, and
J.
Behler
,
J. Chem. Phys.
136
,
064103
(
2012
).
28.
N.
Artrith
and
A. M.
Kolpak
,
Nano Lett.
14
,
2670
(
2014
).
29.
S.
Kondati Natarajan
and
J.
Behler
,
J. Phys. Chem. C
121
,
4368
(
2017
).
30.
Z.
Deng
,
Y.
Mo
, and
S. P.
Ong
,
NPG Asia Mater.
8
,
e254
(
2016
).
31.
A.
Urban
,
D.-H.
Seo
, and
G.
Ceder
,
Npj Comput. Mater.
2
,
16002
(
2016
).
32.
Y. A.
Du
and
N. A. W.
Holzwarth
,
Phys. Rev. B
76
,
174302
(
2007
).
33.
Y.
Deng
,
C.
Eames
,
J.-N.
Chotard
,
F.
Lalère
,
V.
Seznec
,
S.
Emge
,
O.
Pecher
,
C. P.
Grey
,
C.
Masquelier
, and
M. S.
Islam
,
J. Am. Chem. Soc.
137
,
9136
(
2015
).
34.
Y. A.
Du
and
N. A. W.
Holzwarth
,
Phys. Rev. B
81
,
184106
(
2010
).
35.
J.
Lee
,
A.
Urban
,
X.
Li
,
D.
Su
,
G.
Hautier
, and
G.
Ceder
,
Science
343
,
519
(
2014
).
36.
A.
Patil
,
V.
Patil
,
D.
Wook Shin
,
J.-W.
Choi
,
D.-S.
Paik
, and
S.-J.
Yoon
,
Mater. Res. Bull.
43
,
1913
(
2008
).
37.
J.
Kawamura
,
N.
Kuwata
,
K.
Toribami
,
N.
Sata
,
O.
Kamishima
, and
T.
Hattori
,
Solid State Ionics
175
,
273
(
2004
).
38.
N.
Kuwata
,
N.
Iwagami
,
Y.
Tanji
,
Y.
Matsuda
, and
J.
Kawamura
,
J. Electrochem. Soc.
157
,
A521
(
2010
).
39.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
40.
M.
Fuchs
and
M.
Scheffler
,
Comput. Phys. Commun.
119
,
67
(
1999
).
41.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
42.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
43.
H.
Eshet
,
R. Z.
Khaliullin
,
T. D.
Kühne
,
J.
Behler
, and
M.
Parrinello
,
Phys. Rev. B
81
,
184107
(
2010
).
44.
G. C.
Sosso
,
G.
Miceli
,
S.
Caravati
,
J.
Behler
, and
M.
Bernasconi
,
Phys. Rev. B
85
,
174103
(
2012
).
45.
J.
Behler
,
Int. J. Quantum Chem.
115
,
1032
(
2015
).
46.
J.
Behler
,
J. Chem. Phys.
134
,
74106
(
2011
).
47.
S.
Kondati Natarajan
,
T.
Morawietz
, and
J.
Behler
,
Phys. Chem. Chem. Phys.
17
,
8356
(
2015
).
48.
N.
Artrith
,
High-Dimensional Neural Network Potentials for Solids and Surfaces
(
Ruhr-Universität Bochum, Universitätsbibliothek
,
2013
).
49.
G.
Pranami
and
M. H.
Lamm
,
J. Chem. Theory Comput.
11
,
4586
(
2015
).
50.
S.
Sicolo
and
K.
Albe
,
J. Power Sources
331
,
382
(
2016
).
51.
Y. A.
Du
and
N. A. W.
Holzwarth
,
J. Electrochem. Soc.
154
,
A999
(
2007
).
52.
N.
Kuwata
,
X.
Lu
,
T.
Miyazaki
,
Y.
Iwai
,
T.
Tanabe
, and
J.
Kawamura
,
Solid State Ionics
294
,
59
(
2016
).
53.
O. V.
Yakubovich
and
V. S.
Urusov
,
Crystallogr. Rep.
42
,
261
(
1997
).
54.
C.
Keffer
,
A. D.
Mighell
,
F.
Mauer
,
H. E.
Swanson
, and
S.
Block
,
Inorg. Chem.
6
,
119
(
1967
).

Supplementary Material