Predicting highly correlated hydride-ion diffusion in SrTiO₃ crystals based on the fragment kinetic Monte Carlo method with machine-learning potential

The kMC method is described in detail elsewherekMCbook ; kMCbasis01 ; kMCbasis02 . Here, we briefly describe the kMC approach to simulating ionic diffusion. In the jump-diffusion kMC approach, the transition rate ($k_{i}^{j}$) from atom $i$ to atom $j$ can be estimated as

$\displaystyle k_{i}^{j}=A\mathrm{exp}\left[-\frac{E^{\dagger}_{ij}}{RT}\right],$

(1)

where $E^{\dagger}_{ij}$, $R$, $T$, and $A$ are the activation energy, universal gas constant, temperature, and pre-exponential factor, respectively. Summing the respective transition rates gives an estimate of the total rate constant $r^{\mathrm{tot}}$, expressed as

$\displaystyle r^{\mathrm{tot}}=$

$\displaystyle\sum_{i}^{n^{\mathrm{atom}}}\sum_{j}^{n_{i}^{\mathrm{site}}}k_{i}^{j},$

(2)

where $n^{\mathrm{atom}}$ and $n_{i}^{\mathrm{site}}$ are the total number of target atoms in the system and number of nearest neighbor sites for atom $i$, respectively. For example, in the case of an oxide ion in perovskite crystal, $n_{i}^{\mathrm{site}}=8$ for all oxygen sites. The transition probability from $i$ to $j$ ($p_{i}^{j}$) is given by

$\displaystyle p_{i}^{j}=$

$\displaystyle k_{i}^{j}/r^{\mathrm{tot}}.$

(3)

By iteratively selecting the transition from $i$ to $j$ based on the probability $p_{i}^{j}$, the position of target site $i$ is updated. Likewise, the next $n+1$th step simulation time $t_{n+1}$ can be obtained from current time $t_{n}$ and the total rate constant in Eq. 2:

$\displaystyle t_{n+1}=t_{n}-\frac{ln(r^{\prime})}{r^{\mathrm{tot}}},$

(4)

where $r^{\prime}$ is uniform random number from zero to one.

Because the total number of events is $n^{\mathrm{atom}}\times n^{\mathrm{site}}$, the total rate constant $r^{\mathrm{tot}}$ increases cubically with the size of lattice. Therefore, the simulation time step $t_{n+1}-t_{n}$ cubically decreases, and the kMC approach will be limited to relatively small systems. To reduce the computational cost of kMC, several block separation schemes have been developedli2019openkmc ; tada2013 , aiming to achieve good parallel-processing efficiency. We have recently proposed an alternative type of parallelizable scheme involving atom-based fragmentationNAKATA2020109844 , where the transition rate is estimated for each atom $i$ and the maximum rate constant for an atom is defined as

	$\displaystyle r_{i}=$	$\displaystyle\sum_{j\in n_{i}^{\mathrm{site}}}k_{i}^{j},$		(5)
	$\displaystyle r^{\mathrm{max}}=$	$\displaystyle\mathrm{max}\left(r_{i}\right).$		(6)

If we adopt dynamic renormalizationkMCRenormalize , then the transition probability $p_{i}^{j}$ in Eq. 3 can be reformulated as

$\displaystyle p_{i}^{j}=\frac{r^{\mathrm{max}}}{r^{\mathrm{tot}}}\frac{k_{i}^{j}}{r^{\mathrm{max}}}.$

(7)

The transition event selection can then be partitioned into the selection of an atom with $r^{\mathrm{max}}/r^{\mathrm{tot}}$ and the selection of an event for the independent atom $k_{i}^{j}/r^{\mathrm{max}}$. With this renormalization, several ions can be updated simultaneously , and the parallelization of kMC becomes easierNAKATA2020109844 . In this study, we extend this approach by adopting ML potential correction.

The key parameter for kMC simulation is the activation barrier $E^{\dagger}_{ij}$ from $i$ to $j$ defined in Eq. 1. In the case of oxide ion migration in a single SrTiO₃ crystal, the PES is almost flat, as shown by black line in Fig. 1(a), and no correction of the activation barrier of $E^{\dagger}_{ij}$ is required.

However, the flat PES case is quite rare, with the actual PES often being affected by the presence of other types of ions. The main idea in this study is to consider correcting PES changes by using an ML model (the red line in Fig. 1(a)). The transition rate $k_{i}^{j}$ can then be reformulated as

$\displaystyle k_{i}^{j,\prime}=A\mathrm{exp}\left[-\frac{E^{\dagger}_{ij}+\Delta E_{ij}^{\mathrm{ML}}/2}{RT}\right],$

(8)

where $\Delta E_{ij}^{\mathrm{ML}}$ is the potential energy difference between the final and the initial vacancy sites, which can be estimated via the ML model. A similar reformulation of the transition rate can be found in the recent review of Koettgen et al.ceria02 .

As an example, for the case of a vacancy transition in an SrTiO₃ perovskite-type crystal, the vacancy can jump to any of the eight nearest-neighbor sites (see Fig. 1(b),and Fig. 2). First, the transition rates from $i$ to the other sites ($a$, $b$, $\cdots$ $c$) are set from a predefined simulation parameter. The transition rate $k_{i}^{a}$ is then refined using the ML model

$\displaystyle\Delta E_{ij}^{\mathrm{ML}}=f^{\mathrm{ML}}(\mathbf{X}),$

(9)

where $f^{\mathrm{ML}}$ is the ML model and $\mathbf{X}$ is the feature vector used to predict $\Delta E_{ij}^{\mathrm{ML}}$. In this study, the training of the ML model is aimed at developing the structure–energy relationship, with the energy being estimated using first-principles simulation.

The feature vector is a one-dimension array that represents the atomic configuration in the crystal. To show how the vector $X$ is prepared, an example is shown in Fig. 2). In this example, we are trying to construct the feature vector for oxygen vacancy $V_{O}$, depicted in blue. The original position of $V_{O}$ is at the right edge of the periodic boundary condition (PBC) (the yellow line in Fig. 2). We apply the minimum-image convention for $V_{O}$ and the PBC lattice is shifted to the blue line to locate the $V_{O}$ at the center of the PBC (see Fig. 2). The atoms in the oxygen site are then reordered as the first nearest neighbor (red), the second nearest neighbor (green), and the third nearest neighbor (orange). They are labeled down to up and left to right (1, 2, 3, $\cdots$ 10 in Fig. 2). Using this labeling by the minimum image convention, the feature vector $X$ can be filled for any atomic species. In the case of SrTiO₃ perovskite-type crystal, all the oxygen sites are symmetrical, making it possible to use a single ML model $f^{\mathrm{ML}}$ to predict the potential energy difference $\Delta E_{ij}^{\mathrm{ML}}$. By inserting Eq. 9 into Eq. 8, $k_{i}^{j,\prime}$, $r_{i}$, and $r^{\mathrm{max}}$ can be updated and the transition probability $p_{i}^{j}$ can be updated using Eq. 7.

The reference datasets for the ML model were prepared using a standard kMC simulation (i.e., without the ML potential correction) and 1,000 structures were randomly generated in the single kMC simulation run. The system size was 3$\times$3$\times$3 unit cells, as noted in the Computational Details section, and independent kMC simulations were performed for four SrTiO_3-xH_x cases, where $x$ was 0.25, 0.35, 0.45, or 0.60. For each simulation, one or two oxygen atoms were replaced with vacancies, giving eight sets of kMC simulations. Note that the insertion of two oxygen vacancies is at a higher concentration than was used in the actual kMC simulations, but we used both one-vacancy and two-vacancy models for the ML dataset to include the effect of the repulsion potential between vacancies. The size of the dataset for evaluating the structure–energy relation was 8,000 and, using this dataset, we constructed the ML model used to predict the PES in oxyhydrides.

To construct the ML model, the feature vector $X$ in Eq. 9 is prepared using three steps. First, from the trajectory obtained by kMC, the selected transition event (V $\rightarrow$ O, V $\rightarrow$ H) from site $i$ to site $j$ is determined and the atomic labels are assigned using the minimum image convention with respect to the center of the V_O site ( site $i$) (see the Theory and Method section for details). Second, the labels of the atomic species (O, V, H) are replaced by the integer array (0, 1, 2) to obtain the sequential integer vector $X=\left[1,0,0,2,\cdots,2\right]$. To distinguish the transition event from the $i$ to $j$ site, we use the negative integer labels $-1$ and $-2$ for the initial $i$ site and final $j$ site, respectively. Finally, we introduce another integer labeling (using 1 and 2) to distinguish the hydride ion and oxygen ion transitions, with the transition label “1” or “2” being put at the top of the feature vector. The completed feature vector then becomes $X=\left[1,-1,0,-2,2,\cdots,2\right]$, which contains all the necessary information (reaction type, transition site positions, and surrounding atomic species) to determine the change in PES.

The comparison between the ML energy (energy predicted by the ML model) and the density functional theory (DFT) energy (energy estimated by the first-principles simulation) is depicted in Fig. 3. In Fig. 3, the horizontal axis represents the DFT energy and the vertical axis represents the predicted ML energy. The DFT energy denotes the energy difference in the ionic transition from the $i$ site to the $j$ site. As shown in Fig. 3, the ML model reproduces the DFT energy quite well and the root mean square error in Eq. 8 is 0.04 eV, which is less than the error expected for experimentally determined activation barriersliu2019highly . From Fig. 3, we note that the most of the DFT energies are located in a vicinity of zero, indicating that the PES is nearly flat. Furthermore, the comparison between ML and DFT energy shows a similar trend toward positive and negative changes in PES. These results suggest that the ML model developed is adequate for the kinetic simulation of oxyhydrides. We therefore adopt this ML model for predicting the diffusion coefficients of hydride and oxide ions in SrTiO_3-δ.

The hydride ion diffusion coefficients were evaluated for four oxyhydrides (SrTiO_3-xH_x, where x= 0.25, 0.35, 0.45, and 0.60), as noted in the Computational Details section. To evaluate the effect of ML correction on the diffusion coefficients, results for the standard kMC (i.e., without ML) were obtained and these results were compared with those obtained with the help of ML potentials. The results for the diffusion coefficients are shown in Fig. 4, where there is a significant difference between the diffusion coefficients obtained with and without the ML correction. For hydride ion concentrations of $x=0.25$ or $0.35$, the PES accelerates the diffusion of hydride ions significantly, whereas the diffusion coefficient using kMC-ML remains the same for $x=0.45$ and is slightly reduced for $x=0.60$.

We can also note the difference in activation energy with and without PES correction. A summary of the activation barriers with and without ML correction is shown in Table 1. Without the effect of PES correction, the activation barrier monotonically decreased with an increase in the concentration of hydride ions. The largest activation barrier was 0.53 eV for x=0.25, and the smallest was 0.27 eV. The changes of activation barriers can be understood in terms of slow oxygen diffusion. When the hydride ion concentration is small, the diffusion of hydride ions is limited by the surrounding oxygen ions. Therefore, the diffusion of hydride ions is affected by the oxygen ion concentration, with the activation barrier increasing with oxygen ion concentration. However, without PES correction, the estimated activation barriers are slightly larger than those obtained experimentally (See Table 1). Inclusion of ML correction in Eq. 8 changes the activation barrier from 0.42 eV to 0.28 eV for $x=0.35$ and this simulation result with ML correction agrees well with the experimentally derived activation barrier (0.28 eV). Likewise, for $x=0.45$, the simulated activation barrier with ML is 0.32 eV, agreeing well with experiment (0.30 eV), which suggests that the simulation results are reasonable.

Note that the activation barriers estimated via FkMC-ML do not change very much with the variation of hydride ion concentration. However, without ML correction, the activation barrier is significantly changed, as noted above, indicating the importance of the ML potential contribution. For the case of $x=0.25$, the activation barrier of 0.54 eV is close to the oxygen migration barrier in SrTiO₃ (0.6 eV), which suggests that the rate-determining step for hydride ions is oxygen migration. Using ML potentials decreases the barrier to 0.27 eV, which suggest that the inclusion of the PES may affect not only hydride ion migration but also oxygen ion migration, with both hydride and oxygen ion diffusion being accelerated. Because the hydride and oxygen ions interact significantly, the diffusion of oxygen should be analyzed to understand fully how hydride ion diffuse in SrTiO₃ crystals. We therefore also analyzed the kinetics of oxygen diffusion.

To understand fully the kinetics of SrTiO₃ oxyhydride, oxygen ion diffusion was also evaluated for the four hydride ion concentrations. The results for the diffusion coefficients are shown in Fig. 5 and the activation energies are listed in Table 2. Without PES correction, the activation barriers for oxygen ions do not depend on the concentration of hydride ions. The activation barriers are around 0.6 eV, which is the same as for oxygen migration in a single SrTiO₃ crystal.

Inclusion of the ML potential contribution changes significantly the kinetics of oxygen diffusion, particularly for small hydride ion concentrations. When the hydride ion concentration is 0.25, the diffusion coefficient for oxygen using kMC-ML is ten times larger than that obtained without using ML potentials. The contribution of ML potentials to the diffusion coefficient is less noticeable when the simulation temperature is high, which is reasonable because the small changes in potential can be neglected given the increases from thermal fluctuation. Because of the temperature dependence, the apparent activation energy of oxygen ion diffusion decreases from 0.54 eV to 0.35 eV for $x=0.25$. Similar trends can be observed in the simulation results for the other hydride ion concentrations, with the apparent activation energies being decreased by around 0.15 eV using the ML potential contributions for each hydride ion concentration.

To help understand the diffusion of hydride and oxygen ions, the changes in the activation energy for changing concentrations of hydride ions are shown in Fig. 6. In this figure, the difference between the blue and red lines shows the impact of the PES correction, summarizing the impact of interactions between hydride and oxygen ions.

The activation barrier for oxygen ions is uniformly reduced by considering PES changes, with the activation energy being around 0.4 eV (see Fig. 6(b)). By contrast, with hydride ion diffusion, the effect of PES correction does depend on the concentration of hydride ions. There is a reduction in activation energy for small hydride ion concentrations ($x=0.25$ or $0.35$), whereas the activation energy of the hydride ions remains similar for a higher hydride ion concentration ($x=0.45$).

At a low hydride ion concentration ($x=0.25$), hydride ions cannot diffuse without involving the diffusion of oxygen ions. Because the estimated oxygen diffusion is accelerated when considering ML potentials, the estimated diffusion of hydride ions will also be accelerated significantly. For a moderately high hydride ion concentration ($x=0.35$ or $0.45$), hydride ions can diffuse while interacting mainly with other hydride ions. The activation barrier without considering ML potentials will then also become less dependent on the activation barriers for oxygen ions (0.42 eV or 0.32 eV), confirming that slow oxygen diffusion is not the primary rate-determining factor in such cases. Therefore, considering the decreased oxygen activation barrier by using ML potentials does not significantly affect the activation barriers for hydride ions. Likewise, at a high hydride ion concentration ($x=0.60$), a decreased oxygen activation barrier would not accelerate hydride ion diffusion. In contrast to the other cases, using the ML potentials suppresses the diffusion of hydride ions, and the activation barrier is observed to increase from 0.22 eV to 0.27 eV. In summary, the activation barrier for hydride ions is around 0.3 eV, independent of the hydride ion concentration.

In this section, we have investigated the hydride and oxygen ion diffusion mechanisms by analyzing the activation barriers’ dependence on the hydride ion concentration. Considering the inhomogeneous PES created via with ML model, the simulation results concur with the experimental activation barriers obtained by Liu et al.liu2019highly (see Table 1). This offers an insight into why the experimentally determined activation barrier does not depend on the hydride ion concentration. The simulation results indicate that there are two factors determining hydride ion diffusion, namely that increasing the concentration of hydride ions accelerates hydride ion diffusion and that using an inhomogeneous PES accelerates the diffusion of oxygen. These factors affect the rate-determining step, particularly at low hydride ion concentrations.

Finally, the parallel-processing efficiency of FkMC-ML was evaluated for various numbers of CPUs in the range of 64 to 1,000, as noted in the Computational Details section. The results are shown in Fig. 7. In the figure, the red line is the actual computation time, and the black line is the ideally parallelized performance. Ideally, the computational time should decrease as the number of CPUs increases. We obtained a parallel efficiency of 84.92 % using 1,000 CPUs in a comparison with the 108-CPU case. The effective parallel-processing efficiency (the ratio by which the program can be parallelized) was 99.979 %, suggesting that most aspects of the simulation had been parallelized.