Prediction of Born effective charges using neural network to study ion migration under electric fields: applications to crystalline and amorphous Li 3 PO 4

Understanding ionic behaviour under external electric fields is crucial to develop electronic and energy-related devices using ion transport. In this study, we propose a neural network (NN) model to predict the Born effective charges of ions along an axis parallel to an applied electric field from atomic structures. The proposed NN model is applied to Li 3 PO 4 as a prototype. The prediction error of the constructed NN model is 0.0376 e /atom. In combination with an NN interatomic potential, molecular dynamics (MD) simulations are performed under a uniform electric field of 0.1 V/Å, whereby an enhanced mean square displacement of Li along the electric field is obtained, which seems physically reasonable. In addition, the external forces along the direction perpendicular to the electric field, originating from the off-diagonal terms of the Born effective charges, are found to have a nonnegligible effect on Li migration. Finally, additional MD simulations are performed to examine the Li motion in an amorphous structure. The results reveal that Li migration occurs in various areas despite the absence of explicitly introduced defects, which may be attributed to the susceptibility of the Li ions in the local minima to the electric field. We expect that the proposed NN method can be applied to any ionic material, thereby leading to atomic-scale elucidation of ion behaviour under electric fields. IMPACT STATEMENT This study introduces a new computational scheme for analysing ion behaviour in solids under electric fields, through the development of a neural network model to predict the Born effective charges.


I. INTRODUCTION
Ion migration inside various devices, such as all-solidstate batteries and atomic switches [1][2][3], is achieved by applying external forces from applied electric fields.Numerous studies elaborating the stability of materials and the mobility of ions have been conducted using theoretical calculations because these aspects are directly related to device performance.To further advance our understanding of the operating mechanisms of such ionconducting devices, atomic-scale analyses of ionic motion in device operating circumstances, that is, under electric fields, are crucial.
Assuming a linear response, external forces arising from applied electric fields can be estimated simply by multiplying the electric field vector by the valence states of the ions.In electronic state calculations, such as density functional theory (DFT) calculations, the valence states are often evaluated, for instance, using Mulliken charges from the coefficients of atomic orbitals [4] or Bader charges using charge density distributions [5].By contrast, the Born effective charges are defined from the induced polarisation in a periodic system by their atomic displacements (see Fig. 1(a)), or are equivalently defined from the induced atomic forces with respect to the applied electric fields.As our current interest lies in analysing ion motion under applied electric fields, and the latter definition precisely corresponds to the target situation, Born effective charges, rather than static valence states, are the suitable physical quantities to evaluate the external forces acting on the ions.In addition, the Born effective charges can be quantified as the number of each atom without the arbitrariness of the decomposition of the total charges.In most cases, these per-atom quantities are compatible with the computational processes of dynamic calculations using the methods described below.
Recently, atomistic simulations using interatomic potentials constructed via machine learning (ML) techniques have gained increasing attention.Representative methods include the high-dimensional neural network potential (NNP) [6], Gaussian approximation potential [7], moment tensor potential [8], and spectral neighbour analysis potential [9].Numerous studies have demonstrated that ML potentials optimised using DFT calculation data can predict various physical quantities comparable to those of DFT calculations at low computational costs [10][11][12][13].Notably, in their applications to solid electrolyte materials, the predicted ionic conductivities agree well with both the DFT and experimental results [14,15].
To consider the application of these ML potentials for dynamic calculations under applied electric fields, predictive models of ion charges are necessary to evaluate the external forces, as stated earlier.Prior studies have proposed neural network (NN) models to predict the charges of ions [16,17]; however, these models evaluate longrange electrostatic interactions or nonlocal charge transfer through the predicted charge states of ions.Therefore, in this study, we proposed an NN-based model to predict the Born effective charges for given structures.In combination with the conventional NNP, the proposed NN model was applied to dynamic simulations to evaluate the external forces under a uniform electric field.Herein, we employed Li 3 PO 4 as the prototype material, which is commonly used in the research on all-solid-state Li batteries [18][19][20].We verified our scheme of dynamic calculations based on the proposed NN model by evaluating ion behaviour under applied electric fields.

II. METHODOLOGY
Herein, we explain the computational details of the proposed NN model for the Born effective charge predictor.The Born effective charge is defined as follows.
where Ω and e are the cell volume and the elementary charge, respectively; P i and u j are the macroscopic polarisation and atomic coordinates, respectively; F i and E j are the atomic forces and the electric field, respectively; and subscripts i and j represent the x, y, or z directions.In the formalism of NNP, atomic forces can be obtained analytically by applying the chain rule in the following relation.
where U is the total energy and G ν represents the symmetry functions (SFs) [6].Based on the similarities in Eqs.
(1) and (2), the framework of NNP appears to be modifiable as a Born effective charge predictor of a specific direction i (one direction in the 3×3 tensor) by replacing U and F j with − Ω e P i and Z * ij , respectively.Essentially, we confirmed that the above modifications achieved some prediction performance; however, the obtained accuracy was not satisfactory.This inaccuracy can be rationalised by using the scalar quantities from the SFs as the inputs of NN to predict the vector quantities of macroscopic polarisation.
where η and ζ are width parameters; R ij and R ik are the atomic vectors of atom i with j and k, respectively; θ ijk is the angle between atoms i, j, and k at vertex i; α represents either the x, y, or z-coordinates; R s and θ s determine peak positions; and f c is the cutoff function, which is expressed as where R c is the cutoff distance.These functions are invariant to the rotation of the α-axis.In a simplified model, the prediction of Born effective charge are sufficient only for a specific direction along the electric field (e.g., zx, zy, and zz when E z ), as in the following expression (see Fig. 1(b)).
Thus, the proposed NN model requires only a VAF with α = z as its input, which can be achieved by minimal modifications from the original NNP architecture.Note that we may extend the model to predict Born effective charges in the form of tensors in future work.Assuming that the forces vary linearly with the electric field, the external forces acting on the ions can be calculated as follows. ∆F The total forces were considered as the sum of the external forces and values obtained by the conventional NNP.
Thus, simulations of ion dynamics under an electric field can be performed.The loss function of the proposed NN model includes errors in the macroscopic polarisation and Born effective charges in a manner similar to that of the NNP.where N Train and M Train indicate the total amount of data and the total number of atoms, respectively.Here, training was executed by focusing on errors in the Born effective charges (α = 0 and β = 1).We performed atomistic simulations based on the proposed NN model using LAMMPS software [23] with the homemade interfaces.
Next, we describe the training dataset of Li 3 PO 4 used to construct the NNP model.Note that we used parts of the structures and the corresponding total energies and atomic forces of DFT calculations generated in Ref. [14].The dataset includes pristine ( [24] to obtain the Born effective charge tensors.In the DFPT calculations, we used a generalized gradient approximation with the Perdew-Burke-Ernzerhof functional [25], a plane-wave basis set of 500 eV cutoff energy, self-consistent field convergence criterion of 10 All calculations were performed using the Vienna Ab initio Simulation Package [26,27].We used different structural datasets for the two models because the calculated Born effective charges of the strongly distorted structures resulted in unreasonably large values.Additionally, as described above, the proposed NN model predicted three components (one direction of the 3 × 3 tensor) of the Born effective charge for simplicity.Thus, we trained the proposed NN model separately using the Born effective charges in each direction, considering the rotational manipulation of the structures and their tensor values, to enhance the variety of training datasets without performing additional DFPT calculations.

III. RESULTS & DISCUSSION
First, we constructed the NNP using a network architecture of 125 input nodes, two hidden layers with 15 nodes, and one output node, [125-15-15-1], for each elemental species.The root-mean-square errors (RMSEs) of the total energies and atomic forces were 3.34 (2.91) meV/atom and 86.1 (87.9) meV/ Å, respectively, for the randomly chosen 90% (10%) of the training (test) data.The obtained RMSE values were sufficiently small compared with those of other studies using NNP [14,15].Please refer to Fig. S1 for a comparison between the NNP predictions and DFT reference values.The hyperparameters used in the SFs are listed in Tables S1 and  S2.
Next, we constructed the proposed NN model for the Born effective charge predictor.We used the NN architecture of [180-10-10-1], where the RMSEs of the training (randomly chosen 90%) and test (remaining 10%) data were 0.0378 e/atom and 0.0376 e/atom, respectively.Tables S3 and S4 present the hyperparameters used in VAFs. Figure 2 compares the predicted Born effective charges and their DFPT values.Evidently, all the points, including both the diagonal and off-diagonal components, were located near the diagonal lines, thus suggesting fairly good predictions.In addition, the distributions show that the diagonal components of the Born effective charges varied considerably from their formal charges, that is, Li: +1, P: +3, and O: −2, owing to the structural changes.Moreover, the charge states of oxygen underwent the largest variation, despite the fact that Li is a mobile species in Li 3 PO 4 .Furthermore, such variations can be observed in the off-diagonal components of the Born effective charges, although these values were typically small in the crystalline structure: the averages of i̸ =j Z * 2 ij /6 for Li, P, and O are 3.87 × 10 −2 , 2.47 × 10 −3 , and 1.59 × 10 −1 , respectively.In particular, the off-diagonal values of oxygen varied the most between −1 and +1.5.In the following dynamic calculations, we used a uniform electric field of E z = 0.1 V/ Å applied in the z direction.Although the magnitude of the electric field may be excessively large compared with the actual device operating conditions, we chose that value to magnify its effect within a feasible computational time.In addition, the errors of the external forces at the magnitude of this electric field could be estimated to be on the order of meV/ Å, suggesting that the constructed NN model was sufficiently accurate.
Using both the constructed NNP and the proposed NN model, we performed canonical ensemble (N V T ) MD simulations to investigate ion motion under an electric field.For these MD simulations, the temperature and computation time were set to 800 K and 300 ps, respectively.Based on the prior knowledge that Li moves through vacancy sites, we used a crystalline Li 47 P 16 O 64 model, which contained one Li vacancy (V Li ) in the supercell.Note that the size of the simulation model was larger than that of the training dataset, that is, it had a lower V Li concentration.In fact, Li ions seldomly moved in the MD simulations using the pristine model with the above settings or the Li vacancy model at temperatures lower than 800 K.
Figure 3(a) shows the mean square displacement (MSD) of each elemental species, calculated from the trajectories of the MD simulations without an electric field.The MSD of Li increased slightly with MD time, whereas those of P and O remained nearly zero, thus indicating that these elemental species were immobile.In the calculated MSDs under the electric field shown in Fig. 3(b), the MSD of Li increased rapidly, thus indicating Li migration.By contrast, the MSDs of P and O remained small (< 1 Å2 ), as in the case without an electric field.Figures 3(c) and (d) show the MSDs of Li in each direction.We attribute the rapid growth of the total MSD under the electric field to the contributions from the zdirection, which is a physically reasonable result, as it is consistent with the direction of the electric field.In addition, the Li migration paths in crystalline Li 3 PO 4 are not always straight along the z direction, considering that Li moves along the lattice sites via the vacancy hopping mechanism.Hence, the MSDs along the x-and y-directions under an electric field exhibited more fluctuating behaviour compared with the case of without an electric field.
We performed NEB calculations to examine the changes in the potential energy profiles in the presence of electric field.Here, the electrostatic energy that the moving V Li acquires from the electric field, Z * zz E z ∆z, is added to the total energy.For all the NEB calculations, we set the number of intermediate images to six.migration (from images 1 to 8) was reduced from 0.157 to 0.0242 eV by the electric field.By contrast, the potential energy barrier in the opposite direction (images 8 to 1) increased from 0.444 to 0.485 eV.Similarly in the case of Path-2, the potential energy barrier of Li migration from images 1 to 8 (8 to 1) decreased from 0.702 (0.410) to 0.579 (0.402) eV by the presence of electric field.This directionality of the electric field affects the potential energies and facilitates Li migration along the direction of the electric field.In addition, we confirmed that this directionality was almost invisible along a path nearly perpendicular to the electric field (please refer to Fig. S2).The observed decrease in the potential energy barrier and directionality of ion migration agreed with previous NEB calculations of O defects in MgO using the modern theory of polarisation [28].
The atomic coordinates of the migrating Li and the surrounding P and O atoms with or without the electric field in Paths-1 and -2 are shown in Figs.4(d) and (e), respectively.Note that for comparison, we set the initial atomic positions of the migrating Li in the two cases to be identical.We found that the intervals between intermediate images varied according to the presence of an electric field, whereas the paths were similar.Figure 4(f) shows the atomic forces acting on the migrating Li in each NEB image as a function of the z coordinates.Evi-dently, the absolute values of the atomic forces along the z-direction decreased and increased in front of and behind the potential energy barrier, respectively.This indicates a decrease in the barrier height and, subsequently, an acceleration of the Li motion along the z direction.Moreover, we found that the atomic forces along y-direction shifted negatively and positively in the presence of an electric field over Paths-1 and -2, respectively.Because these shifts correspond to the direction of Li motion along y-axis, we considered that the external forces facilitated the movement of Li toward the final positions, whereas their contribution was not as large as those in the z direction.By contrast, the changes in the atomic forces along the x direction were minor because the migrating Li in these paths moved primarily along the z direction, followed by y direction.
Subsequently, we performed additional MD simulations to further examine the effects of external forces along x-and y-directions on Li migration, that is, the contribution from the off-diagonal components of the Born effective charges with respect to the electric field in the z-direction.In these simulations, we considered only the Z * zz , and the Z * zx and Z * zy were set to 0 to exclude their contributions.Thus, we obtained a considerably smaller MSD for Li than that shown in Fig. 3(b) (see Fig. S3).A comparison of the MD trajectory lines shown in Fig. S4 indicates pronounced Li migratory behaviour when all three terms are considered.These results suggest that the off-diagonal components slightly but effectively confined the potential energy surface of the Li migration paths and consequently enhanced Li motion.This enhancement did not appear when the charges were treated as scalar quantities because of the absence of external forces along the x-and y-directions.This also demonstrates the significance of using the Born effective charges to study ion behaviour under electric fields.
Finally, we performed MD simulations with amorphous Li 3 PO 4 under the electric field using the proposed scheme.An amorphous structure was generated using the melt-quench approach, as described in Ref. [14], and the model contained 384 Li, 128 P, and 512 O atoms without including specific defects.Please refer to Fig. S5(a) for the structural image.Here, we set the temperature to 600 K, which is lower than that used in the above cases.Without an electric field, the ions were displaced only to the optimal positions at the early stage of MD, as indicated by the MSDs shown in Fig. S5(b).By contrast, we observed considerably high mobility of Li under an electric field.The obtained MSD value shown in Fig. S5(c) is higher than that shown in Fig. 3(b).Figure 5 shows the trajectory lines of each ion in these two cases.The results clearly show that the ions were immobile without an electric field, whereas the Li ions moved extensively along the z direction in the presence of an electric field.Notably, the P and O ions became relatively mobile in the amorphous structure compared with crystalline Li 3 PO 4 .Furthermore, we found that Li ions moved across the entire region in the amorphous model, whereas Li hopping was restricted to the vacancy sites in the crystal.We also noted that the Li ions moved more pronouncedly at 800 K, as shown in Fig. S6.These results suggest that the Li ions located at the local minima of the metastable amorphous structure susceptibly undergo the electric field effects and readily overcome the mobility barrier.

IV. CONCLUSION
We proposed an NN model to predict the Born effective charges of ions based on their atomic structures.We demonstrated the performance of our proposed NN model using Li 3 PO 4 as a prototype ion-conducting material, where the error in the constructed model reached 0.0376 e/atom.In combination with the conventional NN potential, MD simulations were performed under a uniformly applied electric field.The obtained results indicated an enhanced displacement of Li along the electric field, which is physically reasonable.In addition, we confirmed the lowering of the potential energy barriers from NEB calculations under an electric field.Furthermore, we found that the external forces arising from the offdiagonal terms of the Born effective charges slightly but effectively confined the potential energy surface of the Li migration paths and consequently enhanced Li motion.Finally, we examined the Li behaviour in the amorphous Li 3 PO 4 structure.We found that Li ions located at the local minima were susceptible to electric field effects and readily overcame the mobility barrier.These results suggest that the Born effective charge tensors, depending on the local atomic structures, may be a suitable quantity for a detailed analysis of ion behaviour under external electric fields.

FIG. 1 .
FIG. 1.(a) Schematic of electronic polarisation in Li3PO4 induced by an electric field along z-axis.The yellow (green) clouds depict the increased (decreased) parts of charge density differences, as visualised by VESTA software [29].(b) Schematic of the proposed NN model to predict Born effective charges.

FIG. 2 .
FIG. 2. Comparison between DFPT and NN on the Born effective charges of (a)-(c) training and (d)-(f) test sets.The comparisons are shown separately for each elemental species: (a, d) Li, (b, e) P, and (c, f) O.The light (dark) colours indicate the (off-)diagonal components.

FIG. 3 .
FIG. 3. Calculated MSDs of Li vacancy model (Li47P16O64).The MD simulations with temperature of 800 K (a) without electric field and (b) with Ez = 0.1 V/ Å, where the MSDs are separately shown for each element.The MSDs of Li are separately shown for (c) x and y and (d) z components.

Figure 4 (
a) shows two migration paths of Li, in which Li moved primarily along the z-direction.The potential energy profiles obtained for the two paths are shown in Figs.4(b) and (c).For both paths, we set the potential energy of the initial structure in Path-1 as the reference.In the case of Path-1, the potential energy barrier of Li

FIG. 4 .
FIG. 4. NEB calculations of Li vacancy model (Li47P16O64).(a) Schematic of Li vacancy migrating paths.The potential energy profiles with Ez = 0 and Ez = 0.1 V/ Å, respectively, for (b) Path-1 and (c) Path-2.The atomic coordinates of migrating Li and the neighbouring P and O atoms with Ez = 0 and Ez = 0.1 V/ Å, respectively, for (d) Path-1 and (e) Path-2.(f) Circles, triangles, and squares show the atomic forces of migrating Li for x-, y-, and z-directions, respectively, at each image as a function of its z-coordinate.The open and filled marks show Ez = 0 and Ez = 0.1 V/ Å, respectively.

FIG. 5 .
FIG. 5. Calculated MD trajectory lines of (a, d) Li, (b, e) P, and (c, f) O in the amorphous Li3PO4 model.The upward direction in the figure indicates the z-direction.The MD simulations are performed for 300 ps with the temperature of 600 K (a)-(c) without and (d)-(f) with the electric field.
Li 12 P 4 O 16 ), Li vacancy (Li 11 P 4 O 16 ), and Li 2 O vacancy (Li 22 P 8 O 31 ).For the pristine dataset, we used 14,001 snapshot structures of ab initio molecular dynamics (AIMD) with temperatures of 300, 2000, and 4000 K.The Li vacancy structures comprise 1,656 images from nudged elastic band (NEB) calculations.The Li 2 O vacancy structures comprise 5,000 snapshots of AIMD at 2000 K.In total, 20,657 structures were used to construct NNP.As the training dataset for the proposed NN model, we used 8,000 Li 12 P 4 O 16 snapshots from AIMD calculations at temperatures of 300 and 2000 K, 5,000 Li 11 P 4 O 16 images from NEB calculations, and 4,900 Li 22 P 8 O 31 snapshots from AIMD calculations at 2000 K.For these 17,900 structures, we performed density functional perturbation theory (DFPT) calculations −6eV, and k-point sampling mesh of 6 × 6 × 4 for Li 11 P 4 O 16 and 4×4×2 for Li 12 P 4 O 16 and Li 22 P 8 O 31 .