Accurate meso-scale dynamics by kinetic Monte Carlo simulation via free energy multicanonical sampling: oxygen vacancy diffusion in BaTiO3

ABSTRACT A conceptually accurate method to connect the free energy multicanonical sampling to meso-scale kinetic Monte Carlo (kMC) dynamics is proposed. The required input parameters for kMC simulation are the attempt frequency and activation energy for each event, and the free energy multicanonical sampling enables to obtain the kinetic parameters as a function of temperature, which is the most significant difference from a conventional kMC approach that is based on fixed attempt frequency and activation energy. The present approach is applied to oxygen diffusion in single crystal BaTiO3 including Zn dopant (160 ppm) where an anomaly in the oxygen diffusion is experimentally confirmed; the oxygen diffusion coefficient is slightly dropped at around 1080 K. We carried out 1 μs kMC dynamics in the temperature range of 1020 to 1120 K, and obtained a diffusion anomaly at around 1060 K, which is not obtained in conventional kMC calculations. In addition, the calculated diffusion coefficients using the present approach are in the same order as those of experimental ones, whereas the calculated diffusion coefficients using the conventional method are larger than those of experiment by one order of magnitude at least. The results indicate the advantages of the present approach in comparison with the conventional ones because any assumption and fixation of kinetic parameters are not required in the dynamics simulation.


Introduction
Computer simulations have drawn considerable attention to design new materials from the atomistic point of view. The density functional theory (DFT) [1] is the most successful theory for electronic structure calculations, and combined approaches with molecular dynamics (MD) have revealed atomistic dynamics in materials. Although the DFT-MD approach is one of the most promising computational tools for the investigation of the dynamics without any kinetic information, the approach limits the system size and propagation time in practical computations. Therefore the time acceleration approach by atomistic kinetic Monte Carlo (kMC) [2] simulation has often been adopted to simulate long-term phenomena. For example, the atomistic kMC has been applied to many material simulations such as diffusion in bulk and surface [3][4][5][6][7][8], chemical reaction [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], electrochemical impedance [25], and crystal growth [26,27]. In the rest of this work, we use the term 'kMC' to represent the atomistic kMC.
Despite the successful applications, the system handled in kMC simulation is typically limited to a few thousands of atoms because of several technical problems in reasonable time propagation. Recently, the parallelization based on the domain decomposition with a scheduled execution was applied for the correct time propagation in kMC simulations of much more extended systems, and it enables to simulate the ionic diffusion in a few millions of atoms within standard computational costs [3]. However, we still have another fundamental difficulty in kMC simulations for reliable atomistic dynamics. kMC simulations require a complete set of possible events in a target system. If a possible event is missing, the event never occurs during the simulation, resulting in inaccurate time evolution. Thus, all possible events in a target system should be given in advance of the simulation. Here, the term possible is related to the thermodynamical condition we focus on, and the possible events mean the events which will/might/seldom occur at the target conditions. Now we consider NVT ensemble, where N, V, and T, respectively, mean the number of particles, volume, and temperature of the target system. When we fix N and V, the temperature T is the parameter that makes the meaning of possible clear. Since the target temperature could have a range depending on what we focus on, we have to pick up all possible (i.e. major, minor, and rare) events within a range of target temperatures. Kinetic parameters required for kMC simulations are attempt frequencies and activation energies of all possible events, and the events adopted in kMC simulations are usually limited to what we expect to occur. In addition, the temperature dependence of the kinetic parameters is usually ignored in kMC simulations. However, the reliability of such an approximation is not guaranteed at all. Therefore, accurate kMC simulations require an efficient preprocessing to calculate temperature-dependent kinetic parameters of all possible events within a range of target temperatures.
The extraction of all possible events depending on temperature is itself quite a difficult task. For example, when a system with rugged free-energy landscapes is a target, a lot of computational time would be required to achieve thermal equilibrium. Then, conventional canonical simulations tend to miss many rare but important events for accurate simulations. Recently, some generalized-ensemble methods [28][29][30], such as umbrella sampling [31], parallel tempering (replica exchange) [32][33][34][35][36], adaptive biasing force method [37], metadynamics [38], and the multicanonical (MUCA) ensemble method [39][40][41][42][43][44], have been intensively developed to perform simulations of such a system. Among them, the MUCA ensemble method is especially suitable for the preprocessing of kMC simulation, because the probability distribution function of the MUCA ensemble is a constant within an energy window of interest, which facilitates exhaustive structure samplings within the temperature window. Such an artificial ensemble is made by knowing the density of states of system. In this study, we perform the Monte Carlo simulations with suitable transition rates to find out it. Once we acquire the density of states, we can get the entropy. As a result, we can obtain the statistical quantities in the canonical ensemble, namely free energy, within a range of target temperatures.
The physical/chemical accuracy of MUCA Monte Carlo (MUCA-MC) simulations depends on the energy calculator employed. If we could adopt the firstprinciples electronic structure calculations as the energy calculator, it is absolutely better than the other methods, such as tight-binding approximations and empirical force field approaches. However, from the viewpoint of computational speed, the first-principles calculations are still impractical for the MUCA-MC simulations. Accordingly, an empirical force field approach could be a practical choice in MUCA-MC simulations. Among various force field approaches, we employ the reactive force field (ReaxFF) [45] which can describe various chemical reactions. In general, the parameter set is determined so that the energy from ReaxFF can become equal to that from the first-principles calculations [46][47][48][49]. For the parameter optimization, we have to provide a large amount of structure-related data on crystal formation energy, chemical reaction energy, reaction barrier, bond formation energy, angle potential energy, and torsion potential energy. However, it is not straightforward to prepare such structure-related data, because it is difficult to perform exhaustive structure samplings. Recently, the genetic algorithm has been successfully applied to the parameter optimization for molecular dynamics [50][51][52][53][54]. In this work, we apply the MUCA-MC simulations not only to calculate free energy and kinetic parameters but also to prepare the training data set for the parameter optimization, because the MUCA-MC simulations can generate many structure-related data without any intentional purpose. Using the training data set obtained by MUCA-MC, the ReaxFF parameters are optimized by the random forest regressor.
The route to accurate meso-scale ion dynamics, which we propose in this study, consists of three parts: i) ReaxFF parameter optimizations from the DFT energies, ii) free-energy calculations by the MUCA-MC simulations with the optimized ReaxFF, and iii) ionic dynamics by the parallelized kMC in which each reaction is described by the temperaturedependent kinetic parameters obtained in Step ii). This study is the first report on the multi-scale route from MUCA-MC to kMC simulations.
To demonstrate the powerfulness of the present method, we simulate oxygen vacancy diffusion in meso-scale BaTiO 3 . BaTiO 3 is one of the wellknown compounds as ferroelectric materials which are required in multilayer ceramic capacitors [55][56][57][58][59], ferroelectric random access memories [60,61], and positive temperature coefficient resistors [62,63]. In the electroceramic applications, the oxygen vacancy diffusion is an important property in the design of a robust device because oxygen migration could be an origin of the degradation of the devices. Recently, Kessel et. al. observed the temperature dependence of oxygen vacancy diffusion in BaTiO 3 over a wide high-temperature region [64]. Theoretical investigations have also been reported by classical MD [65,66] and quantum mechanical simulations [67]. However, they mostly focus on the activation barrier of oxygen vacancy diffusion at lower temperatures. Thus, theoretical simulations of oxygen vacancy diffusion at higher temperatures are required, and the temperature dependence of the vacancy migration has to be clarified at the higher temperature region. Therefore, it is worth to employ our meso-scale simulation for ionic dynamics to investigate the temperature dependence of oxygen vacancy diffusion in the low concentration condition. Although the present study is focused on a low concentration of oxygen vacancy for the comparison with experimental results, the present approach can be applied to more complex systems where defects are not diluted, and thus the diffusion is essentially modified by defect-defect interactions.

Optimization of ReaxFF parameter set
The potential energy in the ReaxFF is composed of various contributions: E bond is the contribution of chemical bonds, E lp is the energy penalty related to lone pairs, E over and E under are respectively the energy penalties for over and under coordinations, E val is the contribution related to valence angles (angle energy, penalty energy, and three-body conjugation term), E tors is the contribution related to torsion angles (torsion rotation barriers and four-body conjugation term), E C 2 and E trip are respectively the corrections for C 2 and triple bonds, and E Hbond , E vdW and E Coulomb are respectively the contributions of hydrogen-bond, van der Waals, and Coulomb interactions. (See the supporting information of Ref [49]. for more details). Each contribution has a number of parameters, and thus the global optimization of the ReaxFF parameter set is a very difficult task in general. Recently, such an optimization has been successfully performed by the genetic algorithm [- [50][51][52][53][54], and machine learning (ML) approach [68].
In this study, we adopt the ML approach based on the random forest regressor, as described below. It should be noted that our approach is especially effective when references of the parameter set are available.
As with many optimization problems, the parameter set of ReaxFF is iteratively evaluated until the root-mean-square error (RMSE) of a DFT energy (observed value), E DFT , with respect to a ReaxFF energy (predicted value), E ReaxFF , becomes smaller than a pre-defined threshold. We prepare N p para- to design the regression model for the RMSE on a parameter set and an atomic configuration. The m-th parameter set at k-th iteration P ðkÞ m is given by where 1 and ξ m are respectively the row vectors whose elements are unities and uniform random numbers within [−1,1] for m-th parameter set, δ is the maximum rate of change, and P ðkÀ 1Þ 0 is the optimized parameter set at the previous iteration, more specifically . ., are the parameter subsets for each contribution in Equation (1). Then, the RMSE of m-th parameter set at k-th iteration is written as RMSE ðk;mÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ;I ðkÞ C 1 ; I ðkÞ C 2 ;...;I ðkÞ C Nc ...Þ, which indicates how each parameter is sensitive to the total energy, namely how sensitive it is to the RMSE. Then, we optimize the parameter set through the partial optimizations in order of 'importance', as follows.
For example, let us consider the case that the 'importance' are categorized by three classes: much-important /important/less-important categories. In this case, we start with the optimization for the parameters in the muchimportant category. Assuming that the μ-th parameter set at k-th iteration P ðkÞ μ shows the best RMSE value at that moment and A As a result, we can find from the regression model the parameter set with the best RMSE value among the partially updated parameter sets. After the optimization of the much-important parameters, the regression model is also updated with the partially optimized parameter set. The same procedure is repeated for the parameters in the important and lessimportant categories, and finally we acquire the optimized parameter set at k-th iteration: P ðkÞ 0 .

Multicanonical Monte Carlo simulation
We begin with a short review of canonical Monte Carlo (CA-MC) simulation in atomic system. In the atomic CA-MC simulation, an energy E is determined on a point in the configurational space: Þ is a position of i-th atom and N is the number of atoms, and an atomic position is randomly updated according to the transition rate for a temperature T: where k B is the Boltzmann constant, and E i and E t are the energies of initial and trial configurations in an attempt, X i and X t , respectively. Only an atomic position in X t differs from that in X i by ðδξ x ; δξ y ; δξ z Þ, where ξ x ; ξ y ; ξ z are uniform random numbers within [−1,1] for the displacements in x,y,z directions and δ is the maximum value. The resultant samples satisfy the probability distribution function of the canonical ensemble: where Ω is the density of states of system.
While Ω rapidly grows with E, the canonical weight function e À E=k B T exponentially decays as E increases. As a result, P ca becomes the bell-shaped function with a peak at a specific energy depending on T. In other words, the canonical distribution is very narrow, resulting in a big problem on the simulation. Consider an event that transits from one state to the other one, such as chemical reactions. If the height of free-energy barrier between them is higher than the width of canonical distribution at a certain temperature, the event hardly occurs. Accordingly, the event cannot be captured in time by the CA-MC simulation. Different from the canonical ensemble, the multicanonical (MUCA) ensemble possesses a constant probability within an energy window of interest [39][40][41]. Therefore, if the energy window contains the energy barrier, we could investigate the event precisely by the MUCA simulation. Although the ensemble is artificial, statistical quantities in the canonical ensemble are obtained from the simulation in the last analysis, as described below.
The probability distribution function P muca can be written by where W muca is the MUCA weight function. It should be noted that a variable O is present in the argument of P muca , Ω, and W muca in addition to E. The MUCA ensemble is established as long as there is E in the argument. However, if an appropriate quantity is added to the argument, it facilitates the detailed analysis [69,70]. This extensibility is one of the advantages of MUCA simulations. Actually, we employ a reaction coordinate of an event of interest as O.
There are some rational approaches to obtain the constant probability: P muca ðE; OÞ / const. Here, we employ the approach proposed by Wang and Landau [42,43]. It is immediately found from Equation ( (7)) that P muca becomes a constant value independent of E and O when W muca is inversely proportional to Ω. The MC simulation can be accomplished by the transition rate: Ω is shaped step by step along with the progress of MUCA-MC simulations. In many cases, there are no information on Ω before simulation. Therefore, the development of Ω is the most difficult task in MUCA-MC simulations. For more detail on the calculation procedure, consult the references [42,43].
Once Ω is determined, we can access the statistical quantities in the canonical ensemble. For example, the probability distribution function P and the Helmholtz free energy F along a reaction coordinate are obtained as follows, FðT; OÞ ¼ À k B T ln PðT; OÞ: Here, the subscript 'ca' is dropped for simplicity. According to the harmonic transition state theory, the rate constant of an event that an ion diffuses from A to B site along the reaction coordinate, k A!B , can be evaluated from the landscape of F as follows, where R is the gas constant, F y A!B is the barrier height of the event, and ν A!B is the attempt frequency: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi KðTÞ=m where m is the mass of a diffusing atom and K is the spring constant at the initial state A, which corresponds to the quadratic curvature of F. Using MUCA-MC simulations, we can acquire a rate constant per event which is a key quantity in kMC simulations.

Parallelized kinetic Monte Carlo simulation
Firstly, we briefly introduce the computational algorithm of the conventional kMC simulation for ionic diffusion, in which possible chemical events are ion hopping to neighboring vacancy sites. Considering the case that there are n ion sites with a neighboring vacacy site and the i-th ion diffuses from the ion site, A, to the neighboring vacancy site, B, the kinetics of the ion hopping is characterized by the attempt frequency ν ðiÞ A!B and the free-energy barrier F yðiÞ A!B , resulting in the reaction rate of the hopping k ðiÞ A!B . Then, a total reaction rate R tot is calculated as P n i¼1 k ðiÞ A!B . Using the total reaction rate R tot , we prepare a 'reaction table' whose length is R tot , and the table is divided into small sections numbered from 1 to n in which the individual lengths correspond to the reac- Thus, a reaction with a small/big barrier occupy a long/short length in the table. Based on the reaction table, a reaction to happen at the next time step is determined by a random number in the range from 0 to R tot , and the time step Δt is given as À log γ R tot À � , where γ is another random number in the range from 0 to 1. This is a 'gathering-and-execution' procedure of a single step in kMC. When the number of possible events becomes very huge (i.e. a large size system), the denominator of the time step (i.e. R tot ) could be large, and the single time step will be drastically decreased. This is a problem of the kMC simulation when we extensively extend the system size.
To overcome the difficulty, we have parallelized the kMC algorithm in the domain decomposition scheme where the whole system is divided into several blocks and each block has its own reaction table. Let us consider a case where the whole system is divided into N blocks. The gathering-andexecution procedure is carried out individually in each block at the same time, and we regard the N events happen at the single time step in the parallelized kMC. That is, the total reaction rate R tot is not changed from the non-parallelized execution, but N independent executions are carried out at the same time, which result in N-times faster execution. This is a synchronous approach of the parallelized kMC [3,71,72] in the domain decomposition scheme. A remaining problem of the synchronous approach is an over-selection of events which could happen in a block where the number of possible events is apparently smaller than other blocks. For example, let us consider a system composed of four blocks A, B, C, and D in which blocks A, B, C have 100 events, but D has only one. In a single step of the parallelized kMC, a single event among 100 events is selected in blocks A, B, and C, but only one event in block D is always selected. In a statistical point, the single event selection in block D is the over-selection compared with those in other blocks, and the timing of the single event selection in block D have to be delayed. To avoid the over-selection, we introduce 'Null' events [3,73] which equalize the length of the reaction table of block D to those of blocks A, B, and C by adding no events in block D. When the 'Null' event is selected in block D, nothing happens in the block, which is a reasonable noselection. For a practical use, since the lengths of reaction tables are different with each other, we firstly find a block whose reaction table is the longest one in a system, and individual 'Null' events are introduced in the remaining blocks to realize the same length of the reaction tables in every block.
In addition to the above techniques, we also introduced a scheduled execution for blocks to avoid ion/vacancy vanishing which would occur when an interblock event is selected at both blocks. For example, let us consider the following reaction, , that is, the non-conservation of oxygen happens. Therefore, in this case, the simultaneous executions in blocks A and B have to be avoided by using a scheduled execution. For example, the execution in block B is done after the execution in block A. In this case, the i-th site in block A is changed to a vacancy, and thus the interblock event for block B cannot be a possible event. Along this manner, a scheduled execution for all blocks have to be done in the parallelized kMC in the domain decomposition scheme. The details of the scheduled execution can be found in Ref [3]. The reliability and parallel efficiency of the parallelized kMC are also reported in the same reference for the oxygen ionic diffusion in yttrium stabilized zirconium, a typical oxygen ionic conductor.  Figure 1(a) is the target one to evaluate the temperature-dependent kinetic parameters on the V O diffusion around the Ti atom.

Optimization of ReaxFF parameter set
We extracted the initial Ti-O and Ba-O related parameters from Refs. [74] and [75] respectively, and the Ba-Ti related parameters were tentatively set to the Ti-Ti related parameters in Ref [74]. The ionic charges, which are the parameters in E Coulomb , were fixed to the values extracted from the DFT calculations of the cubic BaTiO 3 . From the above parameters, the initial parameter sets P were prepared by Equation (2) with δ ¼ 0:1. In addition, we prepared the atomic configurations and corresponding DFT energies, E DFT 's. In the structures, 74 structures were generated randomly, 12 structures were designed for oxygen migration, 6 structures for the volume change of BaTiO 3 , and 6 structures for local distortions of Ba-O-Ti. In the partial optimizations in order of 'importance', we set to N q ¼ 14. In the parameter optimizations, scikit-learn code [76] was used for training of random forest regressor, and at the initial stage of the training ReaxFF energies, E ReaxFF 's, were evaluated by the energy calculator implemented in LAMMPS code [77]. It should be noted that the parameter subsets in Equation (3) were limited to those for the four contributions: E vdW , E bond , E val , and E tors , because the four contributions are more responsible for the accuracy than the others and the parameter set used in the previous studies is expected to be not so far from the optimized set for our system. Furthermore, we repeated the optimization until the free-energy barrier of the V O diffusion obtained by the MUCA-MC simulation converges. Therefore, we believe that the updates of only the four contributions are sufficient to obtain the reliable parameter set for our purpose.
The DFT calculations used in the parameter optimizations were performed using the CASTEP code [78,79]. We used the ultrasoft pseudopotentials [80], which can be generated on the fly, to describe the interactions between the ionic core and electrons. The electronic configurations of valence and semicore electrons explicitly treated in the calculations are 5s 2 5p 6 6s 2 for Ba atoms, 3s 2 3p 6 4s 2 3d 2 for Ti atoms, 2s 2 2p 4 for O atoms, and 4s 2 3d 10 for Zn atoms. We employed the exchangecorrelation functional within Perdew-Burke-Ernzerhoftype generalized gradient approximation [81,82]. The cutoff energy for the plane-wave basis set used was taken to be 571.4 eV. We adopted the convergence criterion of 0.03 eV/Å for the geometry optimizations and transition-state searches, and 0.0001 eV for the selfconsistent field calculations of electronic states.

Multicanonical Monte Carlo simulation
We adopted a rotation angle θ of the O atom around the central Ti atom as the reaction coordinate O, as shown in Figure 1(a), which results in exhaustive sampling of the V O diffusion around the Ti atom. In addition to O, we have to specify the maximum displacement δ, and the search domain and its grid points for Ω (in fact, logðΩÞ). Here, δ was set to 0.04 Å, which was determined so that the acceptance/rejection ratio in the MC steps can be almost the same in the CA-MC simulation at 1000 K. The energy window in the search domain was the range from −5.0836 × 10 4 to −5.0240 × 10 4 kJ/mol, which corresponds to the experimental temperature range from 1000 to 1150 K. The angle window was the range from 0 and 90 degrees. The number of grid points was 100 on energy axis and 60 on θ axis. We divided the temperature window into hundredths for the free-energy calculation and thirteenths for the kinetic-parameter calculation.
Our MUCA-MC code has been parallelized with MPI, and 288 CPU cores were used for each MUCA-MC simulation. The computational model is the same as that in the optimization of ReaxFF parameter set, and the energy calculator is ReaxFF implemented in LAMMPS code [77].

Parallelized kinetic Monte Carlo simulation
The computational model for the V O diffusion in single crystal BaTiO 3 consists of 20,479,000 ions (160 × 160 × 160 unit cell of BaTiO 3 ) in which 614 oxygen ions are replaced with V O 's, corresponding to 160 ppm dopant system; as an example, Figure 1(b) shows 30 × 30 × 30 model. The possible dopant elements in BaTiO 3 are Al, Mn, Fe, Zn, Eu, and Y, and we employed Zn as the dopant in this study because it is the major dopant for BaTiO 3 [64]. Ti 4+ ions with the amount of 160 ppm are replaced with Zn 2+ dopant, and thus the concentration of V O is the same as that of the Zn dopant according to the charge neutrality condition. Regarding the temperature dependent kinetic parameters for the V O migration, we executed the MUCA-MC simulation for the V O hopping path shown in Figure 1(a) where the B-site cations are Ti. That is, in the MUCA-MC simulation, we didn't consider the presence of Zn-dopant. In principle, we can execute the same procedure for the V O hopping surrounded by both Ti and Zn, but the amount of the Zn dopant is extremely small and thus we adopted the following approximation: i) the temperature dependences of the kinetic parameters (i.e. hopping barrier and attempt frequency) for the V O hopping are the same between the V O hopping surrounded by Ti only and that by Ti and Zn, and ii) the absolute value of the migration barrier for the V O hopping surrounded by Ti and Zn is determined so as to include the DFT result (0.9635 eV at 0 K).
We performed the kMC simulation using 16 cores in a single node. The temperature range adopted in the kMC simulation was from 1020 to 1120 K. Total number of steps on each kMC simulation was about 1,000,000. For each kMC step, the whole system was divided into 64 sub blocks (4 × 4 × 4), and the event selections were separately performed for each sub block. Thus, the total number of event selections was 64,000,000 times (1,000,000 × 64), which corresponds to about 1 μs. In the calculations of the tracer diffusion coefficient of O, the initial 500,000 steps (about 0.5 μs) were discarded to obtain the value in thermal equilibrium. The tracer diffusion coefficient of the oxygen vacancy D v was calculated every 3,000 steps via the mean square displacement of vacancies, and the tracer diffusion coefficient of oxygen D � is obtained as, where C O and C VO are respectively the concentrations of oxygen and vacancy. It should be noted that the oxygen vacancy concentration adopted in kMC simulation is 160 ppm, and the value is the same as that of the experimental condition.

Accuracy of optimized ReaxFF
In Figure 2, the standardized E DFT or E ReaxFF energies are shown to evaluate the accuracy of obtained ReaxFF parameter. The DFT and ReaxFF data set is standardized as follows, where E i is the E DFT or E ReaxFF energy of geometry i, and μ and σ is the average and standard deviation of E DFT energy. The green filled circle and red open square are the results using initial ReaxFF parameter and optimized ReaxFF parameter, respectively. Black solid line is the ideal line when the E DFT and E ReaxFF energy completely agree well with each other. One can see the significant improvement of ReaxFF parameter  from green to red: After the parameter optimization, the rmse is improved from 1.41 to 0.57, and the Pearson correlation between E DFT and E ReaxFF is also improved from 0.41 to 0.85.
To evaluate the accuracy of optimized ReaxFF parameter furthermore, Figure 3 shows the calculated total energies of E DFT and E ReaxFF in terms of four kinds of structural patterns, along the oxygen migration path (Figure 3(a)), O-Ti-O angle (Figure 3(b)), volume (Figure 3(c)), and Ba-Ti distance (Figure 3(d)). The red open circle are the values of E ReaxFF calculated with the optimized ReaxFF parameter set, and the blue open squares are the values of E DFT . The optimized ReaxFF parameter set reproduce the E DFT energies quite well, and the remaining errors shown in Figure 3 are sufficiently small for practical use; in fact the ranges of errors are almost comparable to other ReaxFF parameter optimizations [83][84][85][86][87][88][89]. Thus, we use the optimized ReaxFF parameter for MUCA-MC simulations.

Kinetic parameters by MUCA-MC simulations
Using the optimized ReaxFF, we performed the MUCA-MC simulations to obtain the density of state (Ω) of BaTiO 3 in which the V O diffusions take place. As a result, the free energy along the reaction coordinate, O ¼ θ, is evaluated from Equation ((10)), leading to the activation barrier and attempt frequency of the event as a function of temperature, F y ðTÞ and v(T). The rate constant necessary for the dynamics of V O by the kMC simulations is calculated from them by Equation (11).
The logarithm of Ω, log(Ω), is shown in Figure 4(a). The horizontal axis is energy (MJ/mol), the vertical axis is θ (degree), and the color bar is log(Ω). One can see that log(Ω) at an energy has a minimum at θ ¼ 45 degrees. This means that the transition state will appear around there, which is clearly observed from the free-energy landscape shown in Figure 4(b). We evaluated the activation barrier, F y ðTÞ, from the landscape. The values are shown in Figure 5(a) and summarized in Table 1. In Figure 5(a), one can see the monotonic increase in the barrier, and the value is 0.6304 eV at 1000 K and 0.6932 eV at 1130 K, respectively. In addition to the activation barrier, we can evaluated the attempt frequency, v(T), from the freeenergy landscape. The result is shown in Figure 5(b). It is found from Figure 4(b) that the curvature around  Figure 1(a), and the color bar is log(Ω). (b) Free energy of V O diffusion in BaTiO 3 crystal as a function of reaction coordinate (angle) and temperature T. The horizontal axis is temperature (K), the vertical axis is the angle (degree), and the color bar is the free energy (kJ/mol).  the bottom of free energy (θ,0 and 90 degrees) becomes larger, resulting in the large spring constant in Equation (12). As a result, the attempt frequency increases with temperature. It should be noted that the temperature dependence can be evaluated by the MUCA-MC simulation. This is one of the most interesting points in our approach. In this case, the temperature dependence is not so large. The difference is only 2%. However, it plays an important role in the dynamics of V O, as shown in next section.

Diffusion coefficient by meso-scale kMC simulations
The calculated results on V O diffusions of BaTiO 3 with Zn dopant (160 ppm) are shown in Figure 6 (a), together with the experimental data; the red circles are the calculated results by using the present method, and the blue squares are the experimental data [64]. Although the slopes of the diffusion coefficients with respect to the temperature show a little difference between the experiment and calculation with the present method, the quantitative differences are less than one order. In fact, a quantitatively good correspondence between them can be confirmed in the linear scale plot (see Figure 6(b)). In addition, a small drop (an anomaly) of the diffusion coefficient at around 1070 K can be observed in both experiment and calculation with the present method, which is originated from a large variation of the diffusion coefficient at around 1060-1070 K, as shown in Figure 6(a) and (c). The large variation of the diffusion coefficients is found only at around 1060 K, and such an anomaly cannot be obtained in the conventional method where the attempt frequency and migration barrier are respectively fixed to be 10 13 /s and 0.67 eV (black circles in Figure 6(a)). That is, an unexpected behavior such as the anomaly of diffusion cannot be obtained by the conventional method, because the slope of the diffusion coefficient is related to the adopted value for the migration barrier in the conventional method. Although the experimental study of BaTiO 3 [64] seems not to take care of the anomaly in their diffusion analysis, our results indicate the anomaly is not an experimental error, and it should be analyzed to figure out the origin of the anomaly, which will be accomplished in a future study. We also confirmed that quantitative difference of the diffusion coefficients by the conventional method from the experimental one is larger than one order. Although the present result based on the combination between the machine learning, multicanonical free energy sampling, and kinetic Monte Carlo is not perfect in terms of diffusion coefficient, we expect the improvements especially in the machine learning will lead to much more satisfactory results, and it will be established as an accurate method.

Conclusion
In this study, we propose a conceptually accurate way to connect the free energy multicanonical sampling to meso-scale kMC dynamics in which atomistic scale information required in kMC is obtained from the firstprinciples simulation without any ambiguity. The required input parameters for kMC simulation are the attempt frequency and activation energy for each event.
In a conventional kMC, the attempt frequency is usually empirically determined to fit the experimental result, or assumed to be a typical number (e.g., 10 −13 /s for ionic conductors), and the activation energy barrier is usually fixed to a single value that is obtained with the firstprinciples calculations. In the present study, on the other hand, we adopted multicanonical free-energy sampling to obtained the attempt frequency and the free energy barrier as a function of temperature to realize an accurate method for large scale ionic diffusion.
As a benchmark target, we investigated oxygen ionic diffusion in a single crystal of Zn-doped BaTiO 3 (160 ppm), because an anomaly in the oxygen diffusion in the system is experimentally confirmed that is clearly deviated from a single slope of diffusion coefficient as a function of temperature. We carried out 1 micro sec kMC dynamics in the temperature range of 1020 to 1120 K, and obtained a diffusion anomaly at around 1060 K, which has not been obtained in the conventional kMC calculations. In addition, the calculated diffusion coefficients using the present approach are in the same order to those of experimental ones; the calculated diffusion coefficients using the conventional method are larger than those of experimental ones by one order at least. The present results suggest the importance of free energy sampling as a function of temperature. Since the present free-energy sampling is executed using the ReaxFF parameters fitted by a machine-learning approach for saving the computational time, the prefect agreement in the diffusion coefficients between the experimental and computational ones is not obtained, but we expect the improvements in the machine learning will lead to satisfactory results for oxygen diffusion, and the present approach will be established as an accurate method for large scale dynamics including chemical reactions.