Phase stability of Au-Li binary systems studied using neural network potential

We firstly performed DFT-based AIMD calculations using the following structures to generate a temporal structural dataset for NNP construction: face centered cubic (FCC) Au, body centered cubic (BCC) Au, BCC Li, FCC Li, ordered FCC Au₃Li (32 atoms/supercell), BCC AuLi (54 atoms/supercell), and rocksalt AuLi₃ (128 atoms/supercell), as well as Au/Li superlattices of BCC Au_0.5Li_0.5 (96 atoms/supercell) and FCC Au_0.5Li_0.5 (160 atoms/supercell). The supercells of these structures are shown in Fig. 1 schematically.

We carried out constant temperature ($NVT$-ensemble) AIMD simulations with 1 fs time step for 1 ps. The initial temperature was set to 300 K and linearly increased up to 2000 K (at 1.7 K/fs) to obtain largely random structural features. Considering the high computational costs of AIMD calculations, further structural generations was conducted using NNP-based MD simulations and subsequent static DFT calculations were performed for the extracted structures.

For DFT calculations, we used the generalized gradient approximation with the Perdew-Burke-Ernzerhof functional PBE , the plane wave basis set (500 eV cutoff energy), and the projector augmented wave method Bloechl-PAW . Brillouin zone integration was performed using the sampling technique of Monkhorst and Pack ($5\times 5\times 5$ and $5\times 5\times 1$ sampling meshes for cubic cells and rectangular cells, respectively) Monkhorst-Pack . We used Vienna Ab initio Simulation Package (VASP) software VASP1 ; VASP2 for all the DFT calculations.

We adopted high-dimensional NNP Behler-PRL-2007 to describe the interatomic potential of Au and Li atoms. In this scheme, the information of atomic structures is encoded as the values of symmetry functions (SFs) and is related to the energy contribution of each atom, $E_{i}$, via the neural network. The total energy of the system is then described by the sum over the energy contributions of the respective atoms ($E=\sum_{i}E_{i}$). Subsequently, the forces along the atomic coordinates $\alpha=x,y,z$ can be expressed as $F_{\alpha_{i}}=-\partial E/\partial\alpha_{i}$ = $-\sum_{\nu}\partial E/\partial G_{\nu}\times\partial G_{\nu}/\partial\alpha_{i}$, where $\nu$ labels the type of SFs. We used the radial ($G_{2}$) and angular ($G_{3}$) SFs in the following forms.

where $\eta$ and $\zeta$ are the width parameters; $\lambda$ and $R_{s}$ determine the peak positions; $R_{ij}$ and $R_{ik}$ are the atomic distances of atom $i$ with $j$ and $k$, respectively; and $\theta_{ijk}$ is the angle consisted of atoms $i$, $j$, and $k$ at the vertex $i$. $f_{\rm c}$ is the cutoff function given by

where $R_{\rm c}$ is the cutoff distance.

The set of hyper-parameters used in the present study are shown in Section S1 and Table S1 in Supplementary Materials. We used 6 and 18 kinds of radial and angular SFs, respectively, for each elemental combination. The combinations for Au atoms include Au-{Au, Li} and Au-{Au-Au, Au-Li, Li-Li} for radial and angular SFs, respectively. Similarly, the combinations of SFs for Li atoms include Li-{Au, Li} and Li-{Au-Au, Au-Li, Li-Li}. Thus, the total number of input nodes (i.e., values of SFs) is $2\times 6+3\times 18=66$ per atom. The SF values were normalized within the range of [$-$1:1] prior to passing through the NN.

Based on the total energies and atomic forces obtained by the DFT calculations, we optimized the weight parameters using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (l-BFGS) algorithm Morales-2011 with the following loss function.

where ${\bf w}$ is the weight parameter vector. $\alpha$ and $\beta$ determine the contribution ratio of energies and atomic forces into the loss function. In the beginning of the NNP training, we evaluated 50 different, randomly assigned weight parameter sets, in which each weight value was chosen within the range of [$-$1:1]. Afterward, we started l-BFGS training with the smallest error weight set using energy differences, i.e., $\alpha=1.0$ and $\beta=0$. Subsequently, we continued the training using both energy and atomic force differences using the setting $\alpha=0.5$ and $\beta=0.5$.

In this section, we will explain the efficient method of NNP construction and structural dataset generation that we developed for the present study. Figure 2 shows the workflow of our NNP construction. First, we constructed a tentative NNP using the structures obtained by the AIMD calculations (please refer to Section II.1 for the computational conditions). Although this tentative NNP is usually less accurate and insufficient for practical use, it was enough to be used for generating additional structures that were added to the training dataset. To obtain these additional structures, we coupled our NNP with Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software LAMMPS1 ; LAMMPS2 using our homemade interface and carried out $NVT$- and $NPT$-ensembles MD calculations for the structures shown in Fig. 1. In the present study, we performed 10-100 ps NNP-MD with a 1 fs time step. The temperatures were set from 300-2100 K at 300 K intervals for the Au case, 300-600 K at 100 K intervals for the Li case, and 300-1200 K at 100 K intervals for the Au-Li alloy cases. The maximum temperatures were chosen based on being several hundreds of Kelvin higher than the materials’ melting temperatures. It is worth noting that rough NNP-MD simulations are often terminated by numerical errors due to structure collapses as mentioned in Ref. Jeong-JPCC-2018 ; Lee-CPC-2019 , wherein the training technique that can construct accurate and numerically stable NNPs was developed.

To select only those structures that have different structural features from the ones already included in the collected dataset, we performed principal component analysis (PCA) based on the SFs of each structure, where $66\times N$-dimensional (where $N$ = number of atoms/structure) local atomic/structural features were reduced to 2-dimensional values, i.e., first and second principal components (PC1 and PC2). Having similar PC values means that those structures have similar structural features. In fact, structures resembling each other can be found frequently in MD trajectories. We calculated SFs for 1% of structures extracted from the NNP-MD trajectories. Then, we chose a specific number of structures by calculating the distances of neighboring points so that the selected PC values were well-scattered (please refer to Section S2 in Supplementary Materials for more details). Figure 3(a-c) shows the PC plots of many of the structures of Au₃Li, AuLi, and AuLi₃, respectively. The plots of the extracted structures are shown in Fig. 3(d-f). The well-scattered plots in all the cases presented suggest that the dataset contained various structures and atomic features (please also refer to Fig. S2 in Supplementary Materials for resulting PC plots). As an example, Figure 4 shows the structures corresponding to the maximum and minimum PC values in rectangular FCC Au_0.5Li_0.5. The maximum PC1 corresponds to the phase separated structure. Conversely, at the minimum PC1 value, the Au and Li atoms were completely mixed. Both the maximum and minimum PC2 show the partial mixture of Au and Li atoms, which is reasonable because their PC1 values were located near the origin. The structures wee more disordered as the PC2 value decreased, although the difference is rather ambiguous compared to PC1 cases. After PCA, we subsequently performed static DFT calculations for the new structures to obtain the total energies and atomic forces. These were then used to update our dataset and retrain our NNP. We performed this cycle iteratively until we could achieve an accurate NNP.

We generated 8696 structures in total, including 1973 Au, 1900 Li, 949 Au₃Li, 992 AuLi, 929 AuLi₃, 971 rectangular BCC Au_0.5Li_0.5, and 982 rectangular FCC Au_0.5Li_0.5 structures. However, we found that at low-Au-concentration alloys, the calculated mixing energies using the NNP seriously deviated from the DFT verifications. Therefore, we further added 589 Au dilute alloy structures, i.e., Au₁Li₃₁, Au₁Li₅₃, and Au₁Li₁₂₇, obtained by the procedure explained above. Finally, our dataset containing 9285 structures covered the entire Au-Li space, suggesting the versatility of the constructed NNP.

We note that Artrith and Behler proposed a method to find structures worth adding to the dataset in the light of structural diversity using NNPs Artrith-PRB-2012 . In their method, NNPs with different network architectures constructed by the same dataset are used to predict energies for a large number of structures (generated for instance by NNP-MD). A large difference in predicted energies of a same structure among NNPs indicates that its atomic features are missing in that dataset. In addition, Li and Ando applied this scheme to the $on$-$the$-$fly$ sampling technique Li-JCP-2019 . Compared with these previous methods, the present approach needed to construct only a single NNP. The lower computational costs for the training process would be the advantage against the above methods.

Using the 9285 structures mentioned in the previous section, we constructed NNPs for the Au-Li binary system. We examined cutoff distances of 5, 7, and 9 Å with the NNs consisting of 2 hidden layers and 10 or 20 nodes per hidden layer. The obtained root-mean-square-error (RMSE) values are summarized in Section S1 and Table S2 in Supplementary Materials. We used a randomly chosen 10% of the structures in the dataset as the test/validation data. We found that all the constructed NNPs predicted nicely both the DFT total energies and atomic forces. The RMSE values are comparable to or even better than the previous studies Behler-PRL-2007 ; Li-JPC-2017 ; Onat-PRB-2018 ; Jeong-JPCC-2018 ; Minamitani-APE-2019 ; Artrith-PRB-2011 . Note that the good prediction performance of the NNP with a cutoff distance of 5 Å suggests the insignificance of long-range interaction in Au-Li systems. Considering the tradeoff between reliability of longer cutoff distances and computational costs, we used the 7 Å cutoff distance and [66-10-10-1] network structure hereafter. The RMSE of total energies and atomic forces in this setting were 1.53 meV/atom and 23.1 meV/Å for the training dataset and 1.46 meV/atom and 22.9 meV/Å for the test dataset, respectively.

Figure 5 shows the comparison between DFT and NNP on the total energies and atomic forces. Both total energies and atomic forces are aligned along the diagonal lines, suggesting that the constructed NNP accurately predicted all structure types considered, i.e., Au, Li, Au₃Li, AuLi, AuLi₃, and dilute-Au alloy. In addition, the learning curves of the total energies and atomic forces are shown in Fig. S1. The mean-square-error (MSE) values monotonically decreased for both total energies and atomic forces, which indicates that no overfitting occured. Note that the preceding 968 ($\alpha=1.0,\beta=0.0$) and 7933 ($\alpha=0.5,\beta=0.5$) iterations were performed ahead of the last training curves shown in Fig. S1.

To verify the accuracy of the constructed NNP, we performed the lattice constant optimization using Au and Li, as well as Au₃Li, AuLi, and AuLi₃ ordered alloys. Figures 6 and S3 show the profiles of total energy as a function of volume per atom obtained by DFT and NNP calculations, together with the fitted Birch-Murnaghan equation of state (EoS). The vertical dashed-dotted and dotted lines correspond to the volume at the energy minima obtained by DFT and NNP, respectively. The DFT and NNP results agreed well on the minimum volume and energy: the maximum differences were smaller than 0.1 ų/atom and 3.78 meV/atom. The bulk moduli were also accurately predicted within a 10% difference, wherein the deviations mainly originated at the compressed structures. Detailed results on fitting to the Birch-Murnaghan equation are given in Section S3 and Table S3 in Supplementary Materials.

Next, we performed phonon calculations for FCC Au, BCC Li, FCC Au₃Li, BCC AuLi, and rocksalt AuLi₃ structures based on DFT and NNP using phonopy software Togo-SM-2015 . We used the same supercell size as shown in Fig. 1 and 0.05 Å displacements to calculate the force constants. Note that we used 0.02 Å displacements for the AuLi₃ case to circumvent the appearance of imaginary frequencies around the X-point.

Figure 7 shows the phonon band structures and the corresponding phonon densities of states (DOS) of FCC Au and BCC Li calculated by DFT and NNP. We found excellent agreement between the DFT- and NNP-obtained phonon modes and the peak structures of phonon DOS. These phonon features could also be found in other studies Lynn-PRB-1973 ; Corso-JPCM-2013 ; Xu-PRL-2014 . Similarly, Figure 8 shows the phonon bands and the corresponding phonon DOS of FCC Au₃Li, BCC AuLi, and rocksalt AuLi₃. The acoustic modes have been well predicted by NNP while the optical modes of Au₃Li show a slight discrepancy around the $\Gamma$-point. However, all the phonon modes match in the AuLi case. This may be attributed to the large number of Au:Li = 1:1 structures included in the training dataset. The phonon bands of AuLi₃ show rather large discrepancies, which is probably due to the relatively sparse sampling of this composition (cf., Fig. S2 especially around the origin). Nonetheless, the phonon dispersion and energy levels were still well reproduced.

It is worth emphasizing that the phonon band structures suggest the high stability of these compounds. Note that we obtained imaginary frequencies for BCC Au and FCC Li phonon bands, which would be reasonable considering the instability of these structures (not shown). This verifies the accuracy of our NNP as well as its capability to reproduce lattice parameters and phonon properties.

As mentioned in the Introduction, the stable Au-Li alloy phases differ among the previous research Pelton-BAPD-1986 ; Bach-EA-2015 ; Bach-CM-2016 ; Yang-JACS-2016 , which arises from the intrinsic complexity of miscible alloy materials and the expensive computational costs. Here, utilizing the constructed NNP, we examined the mixing energies of Au-Li with fine composition grids. We used the ordered structures of Au₃Li, AuLi, AuLi₃, and Au₄Li₁₅ as the base structures (cf., Fig. 1 and Fig. S4 in Supplementary Materials), and Au and Li atoms were randomly replaced for the sake of compositional variations. Note that we replaced Au and Li at intervals of 1 and 4 atoms, respectively, for Au₃Li, AuLi, Au₄Li₁₅, and AuLi₃ base structures. Each structure was first equilibrated by an $NPT$ ensemble MD simulation for 200 ps at $T=$ 500 K. Our choice of $T=$ 500 K is justified by the finding that this value seemed high enough to facilitate the migration of atoms to reach stable configurations. This was then followed by an $NVT$ and subsequently $NVE$ ensembles MD simulations for 200 ps each with the same temperature. Lastly, we fully optimized the structures (including lattice vectors and volume optimization) using the conjugate gradient method with 10^-3 meV and 1 meV/Å of the total energy and atomic forces convergence criteria, respectively. As verification, we checked the total energies of resulting structures using static DFT calculations. The mixing energies were calculated by $E_{\rm mix}=E({\rm Au}_{1-x}{\rm Li}_{x})-\{(1-x)E({\rm Au})+xE({\rm Li})\}$, where $E$ is the energy per atom of the composition designated in the parentheses.

Figure 9 shows the calculated mixing energies using the NNP and static DFT calculations. The squares, circles, triangles, and inverse triangles correspond to Au₃Li, AuLi, AuLi₃, and Au₄Li₁₅ base structures, respectively. The filled symbols indicate the so-called convex hull points, which are connected by the solid lines. The resulting mixing energies show good agreement between NNP and DFT over the whole compositional space. While one exception exists at the convex hull point with $x=0.531$, the extremely small energy difference of the order of 5.80 meV/atom at this point is indistinguishable. From the mixing energy profiles, we found that Au₃Li, AuLi, AuLi₃, and Au₄Li₁₅ consisted of the convex hull, as is the case in the MP. In addition, we found Au_0.469Li_0.531, Au_0.406Li_0.594, and Au_0.344Li_0.656 as the hull points. The structures of these three phases are shown in Fig. 10. The common neighbor analysis method Stukowski-MSMSE-2012 identified that Au_0.469Li_0.531 and Au_0.406Li_0.594 have BCC and FCC crystal structures, respectively. Conversely, Au_0.344Li_0.656 had the laminated structure of two layers of FCC and HCP along the $a$-axis. The radial distribution functions (RDFs) clearly show that Au_0.469Li_0.531 simply replaced a few Au and Li atoms, keeping the AuLi crystal structure. Although the first peaks of Au_0.406Li_0.594 and Au_0.344Li_0.656 were located close to that of Au₃Li, the attenuations of subsequent peaks are indicative of their rather amorphous-like structures. Note that the RDFs are shown in Fig. S5 in Supplementary Materials. In addition to the 7 convex hull points, we also found several points located slightly above the hull over the entire compositional space. These compositions may appear in a specific experimental condition as meta-stable phases.

We also performed the mixing energy calculations for the previously reported stable phases, i.e., Au₂Li, Au₅Li₃, Au₂Li₃, Au₃Li₅, and AuLi₂, starting from the provided space group information Bach-CM-2016 ; Yang-JACS-2016 . The resulting mixing energies of Au₅Li₃, Au₂Li₃, and AuLi₂ were located at 3.38, 6.04, and 4.97 meV/atom below the convex hull of Fig. 9, respectively, while Au₂Li and Au₃Li₅ were 15.1 and 3.28 meV/atom higher than the convex hull, respectively, for the NNP case. Contrastingly, for all five compositions, the mixing energies obtained by DFT were located below the convex hull (please refer to Section S4 and Fig. S6 in Supplementary Materials for details). The mixing energies of Au₂Li, Au₅Li₃, Au₂Li₃, Au₃Li₅, and AuLi₂ were 1.58, 5.88, 4.16, 5.41, and 3.91 meV/atom lower than the convex hull lines at the corresponding composition ratio, respectively. However, the subtle differences in mixing energies in the order of a few meV/atom could be easily buried by circumstances and/or numerical accuracy. Thus, observing a part of the original phase is unavoidable. The discrepancy among previous studies mentioned earlier can be understood from the above reasoning. Note that while considering these phases, Au_0.469Li_0.531 still consisted of the convex hull point.

As demonstrated above, our constructed NNP can be used over the entire range of binary alloy compositions, despite our training dataset including only a limited number of Au-Li ratios, $x$. The precise investigation of mixing energy reveals the compositions that were not present in the previous studies. However, despite the fine intervals of $x$, some crystal structures, inaccessible by the systems used in this work, exist. In fact, the five previously reported crystal structures were not covered in the search of Fig. 9. Further consideration, such as a different number of atoms per supercell, may enable us to draw a more accurate phase diagram. We shall leave this to our future studies.

Next, to clarify the atomic-scale picture of the alloying process, we carried out NNP-MD ($NPT$ at $T=$ 500 K) simulations starting from the phase-separated structure, consisting of 10 atomic layers (32 atoms/layer) of FCC Au and Li using the Au lattice constant (Fig. 11(b)). Figure 11(a) shows the total energy profile along the simulation time, and (b)-(e) show the MD snapshots together with the atomic energy profiles. The colors indicate the magnitude of atomic energies along with the relative atomic positions.

We found that at 4 ps, a single Au atom dissolved into the bulk Li (Fig. 11(c)). This occured due to the slightly larger energy gain of Au dissolving into Li than the Li dissolving into Au, which is corroborated by the mixing energies (cf., Fig. 9). This was followed by continued Au alloying into the Li until around 80 ps, where the total energy reached a plateau (Fig. 11(a)). At this plateau, Li mixed with Au over the entire Li bulk region while 6 layers of Au remained intact (Fig. 11(d)). The cell volume gradually decreased by reducing the length of the $c$-axis, which is attributed to the smaller optimal volume of alloys (cf., Figs. 6 and S3). The alloying was reinitiated at ca. 100 ps with the emergence of a second plateau in the total energy profile. At 200 ps, complete alloying and equilibration (cf., Fig. 11(e)) was reached.

The solution of Au atoms into Li at the initial stage of allying is a probabilistic process, and it proceeds independently at the both Au-Li interfaces. When multiple adjacent Au atoms dissolved into Li, the alloying of the entire interface started from the dissolved region. This corresponds to the steep energy decrease from ca. 4 to 30 ps as shown in Fig. 11(a). In the present case, the dissolution of the lower interface initially took place. The upper interface follows from the inflection point of the energy profile at ca. 12 ps. Note that the initial alloying process using 128 atoms/layer with 10 atomic layers of FCC Au and Li model is shown in Section S5 in Supplementary Materials.

In terms of atomic energies, atoms located at the interface had lower values (Au: $-3.71$ eV, Li: $-2.02$ eV) compared to those of bulk region (Au: $-3.28$ eV, Li: $-1.89$ eV) at the initial structure. Atoms at the second layer from the interface, conversely, had higher values (Au: $-3.23$ eV, Li: $-1.84$ eV). This shows that the Au-Li sharp interface will spontaneously form an alloy. While thermal oscillation caused the atomic energy fluctuation at the bulk region, e.g., 41.0 meV at 4 ps and 51.3 meV at 80 ps, the interface atoms had much lower energies. Once Au atoms mixed with Li atoms, the atomic energies decreased further. A major part of atomic energy change by alloying was absorbed into the Au section. We shall also leave this point for future work.

As demonstrated above, due to the high stability of Au-Li system, Au and Li atoms spontaneously mixed in the simulations within a few hundred picoseconds, which is a noticeably short period compared to the operating time of ionic devices. This alloying property, in other aspects, leads to the stabilization of interfaces of secondary batteries using a Li anode (cf., e.g., Refs. Kato-JPS-2016 ; Kato-SSI-2018 ). On the other hand, controlling the amount of Li introduced into Au must be a significant factor for developing devices using ion conduction, such as VolRAM Sugiyama-aplm-2017 .